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Abstract 



We poorly understand the properties of amorphous systems at small length 
scales, where a continuous elastic description breaks down. This is apparent 
when one considers their vibrational and transport properties, or the way 
forces propagate in these solids. Little is known about the microscopic cause 
of their rigidity. Recently it has been observed numerically that an assembly 
of elastic particles has a critical behavior near the jamming threshold where 
the pressure vanishes. At the transition such a system does not behave as a 
continuous medium at any length scales. When this system is compressed, 
scaling is observed for the elastic moduli, the coordination number, but also 
for the density of vibrational modes. In the present work we derive theo- 
retically these results, and show that they apply to various systems such as 
granular matter and silica, but also to colloidal glasses. In particular we 
show that: (i) these systems present a large excess of vibrational modes at 
low frequency in comparison with normal solids, called the "boson peak" in 
the glass literature. The corresponding modes are very different from plane 
waves, and their frequency is related to the system coordination; (ii) rigidity 
is a non-local property of the packing geometry, characterized by a length 
scale which can be large. For elastic particles this length diverges near the 
jamming transition; (iii) for repulsive systems the shear modulus can be much 
smaller than the bulk modulus. We compute the corresponding scaling laws 
near the jamming threshold. Finally, we discuss the implications of these re- 
sults for the glass transition, the transport, and the geometry of the random 
close packing. 



1. Introduction 



1.1 Anomalous properties of amorphous 
solids 

In the last century, the development of statistical physics revolutionized our 
understanding of matter. It furnished a microscopic explanation of heat, and 
gave a description of different states of matter, such as the liquid or the solid 
state. Later, it explained that sudden transitions between these states can 
occur when a parameter is slowly tuned, despite the microscopic interactions 
staying the same. At the heart of these discoveries lie the concepts of equi- 
librium and entropy. At equilibrium, all the possible states with identical 
energy have an equal probability: this allows to define an entropy, and a 
temperature. Nevertheless, many systems around us are not at equilibrium. 
These can be open systems crossed by fluxes of matter and heat, such as 
biological systems. Another case is glassy systems, such as structural glasses 
or spin glasses, where the characteristic times become so slow that there are 
never equilibrated on experimental time scales. Finally there are also systems 
where particles are too large to be sensitive to temperature, such as granular 
matter. These systems are still poorly understood, and one of the current 
goal of nowadays statistical physics is to explain their original properties, 
and hopefully to find generic methods to describe them. 

We do not have a satisfying description of amorphous systems such as 
structural glasses, colloids, emulsions or granular matter. This is particu- 
larly apparent when one considers the low temperature properties of glasses 
PQ. Their low-temperature specific heat has a nearly-linear temperature de- 
pendence rather than varying as T 3 as would be found in a crystal il]. The 
prevailing explanation for this linear specific heat is in terms of tunneling 
in localized two-level systems [2]: atoms or group of atoms switch between 
two possible configurations by tunneling. This phenomenological model has 
also explain the T 2 dependence of the thermal conductivity at very low tem- 
perature. However, several empirical facts are still challenging the theory 
[31 Hj. Furthermore, after 30 years of research there is yet no accepted pic- 
ture of what these two-levels systems are. At higher temperature, around 



1. Introduction 



7 




Fig. 1.1: Response to a monopole of force in a Lennard- Jones system of 1000 
particles. Black (grey) lines correspond to compressive (tensile) 
stesses. Leonforte et al. [H] 

typically 10 K which corresponds to the Thz frequency range for phonons, 
other universal properties of glasses are not fully understood. In particular 
the thermal conductivity displays a plateau, which suggests that at these 
frequencies phonons are strongly scattered. This effect is significant: for ex- 
ample in silica glass, the thermal conductivity is several orders of magnitude 
smaller than in the crystal of the same composition [Hj. 

Athermal amorphous systems, such as granular matter, also display fas- 
cinating properties, both in their static behavior and in their rheology. The 
following puzzle underlines the subtlety of force propagation in granular mat- 
ter [TU] : the supporting force under a conical heap of poured sand is a mini- 
mum, rather than a maximum, at the center of the pile where it is deepest. 
As we shall discuss in the next Chapter, it has been proposed that in granular 
medium the force propagates differently than in a continuous elastic body 
[TT) IT2] . It turns out experimentally f3J and numerically f5J |H] that an 
elastic-like behavior is recovered at large distances. Curiously enough, the 
cross-over length can be large in comparison with the particle size. Fig. (jl.l|) 
shows the response to a point force in a Lennard- Jones simulations [E] at 
zero temperature. The average response is similar to the one of a continu- 
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ous elastic medium, but near the source the fluctuations are of the order of 
the average. They decay exponentially with distance, with a characteristic 
length of roughly 30 particles sizes. One may ask what determines such a 
distance, below which an amorphous solid behaves as a continuous medium. 
More generally, what length scales characterize these systems? 

The length scales we are discussing might also affect the rheology of gran- 
ular matter. An interesting question is how grains flows, or how they compact 
For example if a layer of sand is inclined, an avalanche is triggered. In- 
terestingly the angle 9 of avalanche appears to be controlled by the width h 
of the granular layer. 9 decreases when h grows when h is smaller than of 
the order of ten particle sizes. Similar length scales also appear in the spatial 
correlations of the velocities of grains in dense flows [T7j . 

A particularity of the amorphous state is that it is not at equilibrium. 
Consequently the properties and the microscopic structure of these systems 
depend much on their history. For example if a granular pile is made by a 
uniform deposition, rather than by pouring sand from the top, the supporting 
force does not display the minimum discussed above at the center of the pile, 
but rather a flat maximum. Often amorphous solids are obtained from a fluid 
phase by varying some parameters such as temperature, density or applied 
shear stress until the system stops flowing: this is the jamming transition 
|THj . As the dynamics greatly slows down once this transition is passed, the 
structure of amorphous solids does not to differ too much from the marginally 
stable state at the transition. Thus a better understanding of the microscopic 
features of amorphous solids requires a better knowledge of the jamming 
mechanisms. It is a hard and much studied problem. When a glass is cooled 
rapidly enough to avoid crystallization, the relaxation times rapidly grow. In 
some cases the relaxation times follow an Arrhenius law with temperature; 
such glasses are called "strong". If the relaxation times grow faster, the 
glass is "fragile" ^H] • There is no available theory to compute quantitatively 
the temperature dependence of the relaxation times, and to decide a priori 
which glasses are strong or fragile. Recently it was observed numerically and 
experimentally that the relaxation in the super-cooled is very heterogeneous 
and involves rearrangements of particles clusters j^BI^D]- Althought several 
models of the glass transition predict such heterogeneities, see e.g. |22] and 
references therein, their cause and nature is still a much debated question. 

Although they can lead to collective dynamics, most of the spatial models 
of the relaxation near the jamming threshold have purely local rules. This is 
the case for example for kinetically constrained models where particles 
are allowed to move individually if their direct neighborhood satisfies some 
specific conditions. The starting point of the present work is the following 
remark: the stability against individual particle displacement is much less 
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demanding than the stability toward collective motions of particles. In d 
dimensions, d+1 neighbors are sufficient to pin one particle. As we shall dis- 
cuss in details later, Maxwell showed that 2d contacts per particle on average 
are necessary to guarantee the stability of a solid [2A j . The fact that the cri- 
terion of stability is non-local suggests that it is so for the minimal motions 
responsible for the relaxation. In any case, this underlines the importance 
of understanding what guarantees the stability of an assembly of interacting 
particles. In what follows we aim to furnish a microscopic description of the 
rigidity characterizing amorphous solids. 

The informations about the rigidity of a solid against collective particle 
motions are contained in the density of the vibrational modes D(u): a system 
is stable if there are no unstable modes. In a continuous isotropic elastic 
medium, the invariance by translation implies that the vibrational modes are 
plane waves. As a consequence in three-dimensions the density of vibrational 
modes D(uj) follows the Debye law D(u) oc uj 2 . By contrast, at low frequency 
all glasses present an excess of vibrational modes in comparison with the 
Debye behavior. This excess of vibrational modes is the so-called "boson 
peak" 1 which appears as a maximum in D{uj)/uo 2 . It is observed in particular 
in scattering experiments, see Fig, (jl. 2)1 . The frequency of the peak lies in 
the terahertz range, that is typically between cj^/IO and cj^/IOO, where uoo 
is the Debye frequency. 

Several empirical facts suggest that the presence of these excess modes is 
related to many of the original properties of amorphous solids. In most glasses 
|23 ESI EDI EOj, with some exceptions as silica [SI], this boson peak shifts 
toward zero frequency when the glass is heated, as shown in Fig (jl.2j) . Even- 
tually the peak reaches zero frequency, as it has been observed numerically 
jS2] and empirically [2*5] . This suggests that in some glasses the correspond- 
ing modes take part in the relaxation of the system HH] • The presence of 
the boson peak also affects the low-temperature properties of glasses. The 
plateau in the thermal conductivity appears at temperatures that correspond 
to the boson peak frequency [T] , which suggests that the excess-modes do not 
contribute well to transport. Furthermore, since these modes are soft, one 
expects that their non-linearities are important. Hence these modes may 
form two- levels systems IHfi] . Finally, as the linear response to any force 
or deformation can be expressed in terms of the vibrational modes, it is rea- 
sonable to think that the boson peak affects force propagation. The recent 
simulations from which Fig. (jl.l|) is taken show that the length scale that 
appears in the response to a point force also appears in the normal mode 

1 The term "boson peak" was introduced because the amplitude of the scattering peak 
varies according to the Bose-Einstein factor at low temperature. 
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analysis: only for larger system sizes the lowest frequency modes are the one 
expected from a continuous elastic description 

This excess of vibrational modes has been studied with diverse approaches. 
There are phenomenological models also dealing with two-levels systems such 
as the "soft potential theories" EH EE] , which assumes the presence of 
strongly anharmonic localized soft potential with randomly distributed pa- 
rameters. A second approach consists in studying the vibrations of elastic 
network with disorder. Simulations of a harmonic lattice with a random dis- 
tribution of force constants (SHI HOj exhibit a density of states qualitatively 
similar to what is observed with glasses in scattering experiments. From the 
theoretical point of view, models that assume spatially fluctuating elastic 
constants [JT] show an excess of modes whose frequency decreases with the 
amplitude of the disorder. Recently further developments were proposed us- 
ing the euclidean random matrix theory HH1 HI], where an assembly of 
particles at infinite temperature is considered. The density of states corre- 
sponds to the spectrum of a disordered matrix, the dynamical matrix. When 
the density p is infinite, the system behaves as a continuous medium. To 
approximate the density of states at finite p one uses perturbation theory in 
the inverse density. This leads to an excess of modes whose frequency goes 
to zero as the density decreases toward a finite threshold, and furnishes sev- 
eral exponents that describes the density of states at this transition. A third 
approach uses the mode coupling theory (MCT). MCT models the dynamics 
of supercooled liquids, and predicts a glass transition at finite temperature 
where the relaxation time diverges. In the glass region, the structure does not 
fully relax for any waiting time. Nevertheless if the mode coupling equations 
are used into this glass region, the dynamic that appears is non-trivial. It can 
be interpreted in terms of harmonic vibrations [45, around a frozen amor- 
phous structure. The corresponding spectrum displays an excess of modes 
that converges toward zero frequency at the transition. 

These different approaches describe the presence of an elastic instability, 
and make predictions on how the density of states behaves with some pa- 
rameters when this instability is approached. Nevertheless they have several 
drawbacks. In these models, disorder is the main cause for the excess den- 
sity of states. This is inconsistent with scattering data that show that some 
crystals also display excess vibrational modes EH SHI EH] ■ In particular 
silica, the glass with one of the strongest boson peak, as a density of states 
extremely similar to the crystals of identical composition and similar density, 
see Fig. (|8.2jh as we shall discuss in more details in Chapter |H1 Thus the cause, 
and more importantly the nature of these excess modes are still unclear. In 
particular one may ask what microscopic features determine the vibrational 
properties of amorphous solids at these intermediate frequencies, and what 
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Energy (meV) 



Fig. 1.2: Reduced density of states (g(E) ~ D(u) following our notation) of 
collective motions in toluene, ethylbenzene, dibutylphthalate, and 
glycerol glasses. Arrows indicate the energy of the boson peak 
estimated from the data at lowest temperature. Chumakov et al., 

EDI- 
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are the signatures of marginal stability at a microscopic level. In what follows 
we attempt to answer these questions in weakly-connected amorphous solids 
such as an assembly of repulsive, short-range particles. Then we argue that 
this description also applies to other systems, such as silica glass or parti- 
cles with friction. Finally we derive some properties of colloidal glasses and 
discuss the possible implications of this approach for the glass transition. 

1.2 Critical behavior at the jamming 
transition 

Recently, C.S O'Hern, L.E Silbert et al. EI] exhibited a system whose 
vibrational properties are dramatically different from a conventional solid. 
They simulate frictionless repulsive particles with short range interactions 
at zero temperature. The authors consider soft spheres. For inter-particle 
distance r < a, the particles are in contact and interact with a potential: 



where a is the particle diameter and e a characteristic energy. For r > a the 
potential vanishes and particles do not interact. Henceforth we express all 
distances in units of a, all energies in units of e, and all masses in units of the 
particle mass, m. The simulations were done for a = 2 (harmonic), a = 5/2 
(Hertzian contacts) and a = 3/2. When the packing fraction is low, such 
system is in a gas phase, and the pressure p is zero. At high packing fraction, 
it forms a solid and has a positive pressure. There is a transition between this 
two phases where the pressure vanishes: this is the jamming transition. At 
that point the density of states behaves as a constant instead of the quadratic 
dependence expected for normal solids, see Fig fjOj) . The authors also study 
the solid phase when the pressure decreases toward zero. They find that 
the jamming transition acts as a critical point: the microscopic structure, 
the vibrational modes and the macroscopic elastic properties display scaling 
behaviors with the pressure p or with <fi — <f) c , where <p c is the packing fraction 
at the transition. In three dimensions <p c — *> 0.64 which corresponds to the 
random close packing 2 . Concerning the structure, the coordination number 
z, which is the average number of contacts per particles, is found to follow: 



2 The parameter tfi — tfi c is somewhat less natural than the pressure because (f) c can vary 
from sample to sample. The distribution of <j) c converges to a well-defined value only when 
the number of particle N diverges. Nevertheless, the parameter cf) — <j) c has the advantage 
of being purely geometrical, and following 50 we should use it in most cases. 
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independently of the potential, where z c = 2d, and d is the spatial dimen- 
sion. This singular increase of the coordination was already noticed in [52] . 
Another striking observation is the presence of a singularity in the pair cor- 
relation function g(r) at the jamming threshold. g(r) has an expected delta 
function of weight z c at a distance 1 that represents all particles in contact. 
But it also displays the following singularity: 

9ir) ~ -jL= (1.3) 

which indicates that there are many particles almost touching. Again this is 
independent of the potential. This property was observed in other situations 
|53j . Note that Eq. (jl.2j) and (jl.Hj) are related. A small affine compression of 
the configuration is equivalent to an increase of the particle diameter of an 
amount ~ ^r^ 2 - At the jamming threshold this would lead to an increase in 
the coordination number: 



1 1 <t>— <S>C 



i 



z — z c ^4:TT c g(r)r dr ~ / dr ~ (0 — C ) 2 (1.4) 

Ji Ji vr — 1 

as observed. 

As we mentioned, these simulations also reveal unexpected features in the 
density of states, D(u): (a) As shown in Fig (jOj) . when the system is most 
fragile, at p — > 0, D(u) has a plateau extending down to zero frequency with 
no sign of the standard u 2 density of states normally expected for a three- 
dimensional solid, (b) As shown in the inset to that figure, as p increases, 
the plateau erodes progressively at frequencies below a frequency u*, which 
scales in the harmonic case as: 

u* ~ 5z (1.5) 

(c) The value of D(uS) in the plateau is unaffected by this compression, (d) 
At frequency much lower than u*, D(oS) still increases much faster with u 
than the quadratic Debye dependence. 

Finally, it was conjectured by Alexander that the elastic moduli 
should scale at the jamming transition. Such scaling properties where ob- 
served in emulsions near the jamming transition but the presence of 
noise in the measure of the packing fraction makes the exponent hard to 
identify. In jSO], the bulk modulus B, the shear modulus G and the pressure 
are found to follow: 

(1.6) 
(1.7) 
(1.8) 



B ~ ( 




G ~ ( 




p ~ ( 
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Fig. 1.3: D(uS) vs angular frequency uj for the simulation of Ref [SD]- 1024 
spheres interacting with repulsive harmonic potentials were com- 
pressed in a periodic cubic box to volume fraction <p, slightly above 
the jamming threshold <fi c . Then the energy for arbitrary small 
displacements was calculated and the dynamical matrix inferred. 
The curve labeled a is at a relative volume fraction — C = 0.1. 
Proceeding to the left the curves have relative volume fractions 
10~ 2 , 10~ 3 , 10" 4 , 10" 8 , respectively. Inset: Scaling of uj* vs 8z. uj* 
for each (0 — <f> c ) is determined from the data in the main panel 
as the frequency where D(u) is half of the plateau value. 5z vs 
(<f) — <p c ) is obtained from the scaling measured in jHO]- The line 
has slope 1. 
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These results raise many questions. Among others: (i) the vibrations of 
a normal solid are plane waves; what are the vibrations of a random close 
packing of elastic spheres? (ii) How is the behavior of the structure and 
the density of states related? This critical behavior is a stringent test to 
such theories, (iii) How does the microscopic structure, for example the 
coordination number, depend on the the system history? (iv) What are the 
different elastic properties of this system, that behaves almost as a liquid at 
the transition, as the shear modulus becomes negligible compare to the bulk 
modulus? For example, how does it react to a local perturbation? 

1.3 Organization of the thesis 

At the center of our argument lies the concept of soft modes, or floppy modes. 
These are collective modes that conserve the distance at first order between 
any particles in contact. They have been discussed in relation to various 
weakly-connected networks such as covalent glasses [221 E0]> Alexander's 
models of soft solids jHE], models of static forces in granular packsjHH ED) 
and rigidity percolation models, see e.g. [01]. As we will discuss below, they 
are present when a system is not enough connected. As a consequence, as 
Maxwell showed |24j . a system with a low average coordination number z 
has some soft modes and therefore is not rigid. There is a threshold value 
z c where a system can become stable, such a state is called isostatic. As 
we shall discuss later, this is the case at the jamming transition, if rattlers 
(particles with no contacts) are excluded. There are no zero-frequency modes 
except for the trivial translation modes of the system as a whole. However, 
if any contact were to be removed, there would appear one soft mode with 
zero frequency. Using this idea we will show in what follows that isostatic 
states have a constant density of states in any dimensions. When z > z c , the 
system still behaves as an isostatic medium at short length scale, which leads 
to the persistence of a plateau in the density of states at high frequency. 

The second concept we use is at the heart of the work of Alexander on 
soft solids [HO] • In continuum elasticity the expansion of the energy for small 
displacements contains a term proportional to the applied stress (that we 
shall also call initial stress term following [HO]), as we shall discuss in the 
next Chapter. It is responsible for the vibrations of strings and drumheads, 
but also for inelastic instability such as the buckling of thin rods. Alexander 
pointed out that this term has also strong effects at a microscopic level in 
weakly-connected solids. For example, it confers rigidity to gels, even though 
these do not satisfy the Maxwell criterion for rigidity. We will show that 
althought this term does not affect much the plane waves, it strongly affects 
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the soft modes. In a repulsive system of spherical particles it lowers their 
frequency. We shall argue that this can change dramatically the density of 
states at low frequency, as it will be confirmed by a comparison of simulations 
where the force in any contact is present, or set to zero. We show that these 
considerations lead to a inequality between the excess connectivity 5z = z—z c 
and the pressure that guarantees the rigidity of such amorphous solids. This 
relation between stability and structure will enable us to discuss how the 
history affects the microscopic structure of the system. In particlar, we shall 
argue that the preparation of the system used in [SOJ leads to a marginally 
stable state, even when > <f> c . This will account for both the scaling of the 
coordination, and for the divergence of the first peak in g(r) at the random 
close packing. 

A distinct and surprising property of the system approaching the jam- 
ming threshold is the nature of the quasi-plane waves that appear at lower 
frequency than the excess of modes. The peculiar nature of the transverse 
waves already appears at zero wave vector, as the shear modulus becomes 
negligible compared to the bulk modulus near jamming. If the response to 
a shear stress were be a perfect affine displacement of the system, the cor- 
responding energy would be of the same order of the energy induced by a 
compression. As this is not the case, this indicates the presence of strong 
non-affine displacements in the transverse plane waves. To study this prob- 
lem we shall introduce a formalism that writes the responses of the system 
in terms of the force fields that balance the force on every particle. This 
enables us to derive the scaling of the elastic moduli, and to compute the 
response to a local perturbation at the jamming threshold. We show that 
this response extends in the whole system. 

We study how these ideas apply to real physical systems, such as granular 
matter, glasses and dense colloidal suspensions. In granular matter friction 
is always present. We derive the equation of the soft modes with friction. 
The main difference with frictionless particles is that rotational degrees of 
freedom of grains now matter, but our results on the vibrational modes and 
on the elastic properties are unchanged. Then we discuss the case of glasses. 
In these systems the coordination number is not well defined, as there are 
long range interactions such as Van der Waals forces. We show that if the 
hierarchy of interactions strengths is large enough, our description of the 
boson peak still applies. In particular we argue that the boson peak of silica 
glass corresponds to the slow modes that appear in weakly connected systems. 
Our argument also rationalizes why the crystal of identical composition and 
similar density, the crystobalite, has a similar density of state. We propose 
testable predictions to check if the same description holds for Lennard- Jones 
systems. Finally we study dense hard sphere liquids. Our main achievement 



1. Introduction 



17 



is to derive an effective potential that describes the hard sphere interaction 
when the fast temporal fluctuations are averaged out. Our effective potential 
is exact at the jamming threshold. It allows to define normal modes and to 
derive several properties of a hard sphere glass. In particular it implies that 
the jamming threshold act as a critical point both in the liquid and in the 
solid phases. This suggests original relaxation processes. 

The thesis is organized as follows. Chapter El is introductive: we define 
rigidity and soft modes and discuss how these concept were used to study 
covalent glasses, force propagation or gels. In the Chapter 0] we use a simple 
geometric variational argument based on the soft modes to show that isostatic 
states have a constant density of states. The argument elucidates the nature 
of these excess-modes. In the ChapterHJwe compute the density of states when 
the coordination of the system z = z c +5z increases with the packing fraction. 
At that point we neglect the effect of the applied stress on the vibrations. 
This approximation corresponds to a real physical system: a network of 
relaxed springs. We show that such system behaves as an isostatic state for 
length scales smaller than I* ~ bz~ l . This leads to a plateau in the density 
of states for frequency higher than uo* ~ 8z. At lower frequency, we expect 
that the system behaves as a continuous medium with a Debye behavior, 
which is consistent with our simulations. We extend this result to the case of 
tetrahedral networks. In ChapterElwe study the effect of the applied pressure 
on D{uj). We show that althought it does not affect much the plane waves, 
the applied stress lowers the frequency of the anomalous modes. We give a 
simple scaling argument to evaluate this effect, and we discuss its implication 
for the density of states. Incidentally this also furnishes an inequality between 
5z and the pressure which generalizes the Maxwell criterion for rigidity. We 
discuss the different length scales that appear in the problem. In Chapter 
El we discuss the influence of the cooling rate and the temperature history 
on the spatial structure and the density of states of the system. We show 
that the scaling of the coordination and the divergence in g(r) are related to 
the marginal stability of the system of [SU]. Some elastic properties of this 
tenuous system are computed in Chapter [7|, in particular the elastic moduli. 
Chapter |S1 is devoted to the applications of our arguments to granular matter 
and glasses. This approach explains the qualitative shape of the density 
of state of silica. In Chapter 03 we derive an effective potential for hard 
spheres, and compute some properties of a hard spheres liquid near the glass 
transition. To conclude in chapter 11 we discuss the possible applications of 
these ideas to the low temperature properties of glasses, the glass transition 
and the rheology of granular matter. 



2. Soft Modes and applications 



2.1 Rigidity and soft modes 

More than one century ago, Maxwell [21], working on the stability of en- 
gineering structures, studied the necessary conditions for the rigidity of an 
assembly of interacting objects. His response is as follows: consider for ex- 
ample a network of N point particles connected with iV c relaxed springs of 
stiffness unity in a space of dimension d. The expansion of the energy E is: 

5E = W\m 3 -5R i ).n l] f (2.1) 

where the sum is taken on every couple of particles in contact (ij), is 
the unit vector going from i to j, and 5Ri is the displacement of particle i. 
It is convenient to express Eq. (|2.1|) in matrix form, by defining the set of 
displacements SRi...SRn as a diV-component vector |5R). Then Eq. (j2.1|) 
can be written in the form: 

5E=(5R\M\5K) (2.2) 

The corresponding matrix Ai is known as the dynamical matrix , see foot- 
note 3 for an explicit tensorial notation. The 3N eigenvectors of the dynamical 
matrix are the normal modes of the particle system, and its eigenvalues are 
the squared angular frequencies of these modes. A system is rigid if it has 
no soft mode, which are the modes with zero energy. Since Eq. (|2.1j) is a sum 
of positive terms, such modes satisfy: 

(SRi — 5Rj).fiij = for all N c contacts (ij) (2.3) 

This linear equation defines the vector space of displacement fields that con- 
serve the distances at first order between particles in contact. The particles 
can yield without restoring force if their displacements lie in this vector space. 
Fig ()2.1|) furnishes an example of such mode. Note that Eq. (|2.3|) is purely ge- 
ometrical and does not depend on the interaction potential. 

3 M. can be written as an N by N matrix whose elements are themselves tensors of rank 
d, the spatial dimension JAij = — f ® ^ij + S</> ™U ® n>U where 8uj\ = 1 

when i and j are in contact, the sum is taken on all the contacts I with i. 
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Fig. 2.1: Illustration of a rigid and a floppy network made of 4 particles 
interacting with springs. Nd — d(d+l)/2 = 5 contacts are required 
to rigidify the structure. The red arrow of the floppy network 
indicates the soft mode that appears when one contact from the 
rigid system is removed. 

Maxwell noticed that Eq. ()2.3|) has N c constraints and Nd — d(d + l)/2 
degrees of freedom if we substract global translations and rotations. Each 
equation restricts the dN — d(d + l)/2-dimensional space of |5R) by one 
dimension. In general, these dimensions are independent, so that the number 
of independent soft modes is dN — d(d + l)/2 — N c . A rigid system must 
not have any soft modes, and therefore as at least as many constraints as 
it has degrees of freedom. For a large system this yields for the average 
coordination z = 2N C /N: 

z>2d (2.4) 

This is the Maxwell criterion for rigidity. It is important to note that it is 
a global criterion, as it discusses the stability toward collective motions of 
particles. A local criterion that treat the motions of single particles only 
leads to z > d + 1. As we shall see, some systems live on the bound of 
Eq. (j2.4|) . they are called isostatic. 

If this criterion has been known for quite a long time, it is only in the 
last decades that the concept of soft mode of non-rigid systems was used to 
study gels, covalent glasses and force propagation in granular matter. In the 
present Chapter we discuss these ideas. 

2.2 Soft modes and force propagation 

Recently several theories were proposed to describe the response to forces 
in sand. Some E] have argued that granular matter requires a new 
constitutive law: they postulate a linear relation between the components of 
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the stress tensor, the "null stress law". This leads to a continuum mechanics 
different from elasticity, where forces propagate along favored directions. In 
experiments, this description breaks down on large length scales Ej 
where an elastic-like behavior is recovered. Nevertheless, as we shall see 
now, this theory was further justify for frictionless grains |oTl 13HI EH] using 
the concept of soft modes. 

There is an obvious connection between soft modes and forces: soft modes 
do not have restoring force. Therefore a non-rigid system can resist to an 
external force |F) = {Fi}, where % labels the particles, only if the external 
force is orthogonal to each soft mode (3, that is: 

(F\6R P ) = £ F ■ 5R^ = (2.5) 

i 

If Eq. ()2.5|) is not satisfied, the system yields along the soft modes, which are 
thus the directions of fragility of the system. 

To apply this idea to granular matter, the starting point is the follow- 
ing remark: an assembly of frictionless hard spheres (or equivalently elastic 
spheres at the jamming transition where the pressure vanishes) is exactly 
isostatic, as was shown in particular in [23 E3 EH] and confirmed in the 
simulation of |5()j . The argument for hard spheres is as follows: on the one 
hand, it is a rigid system and therefore must satisfy the bound of (|2.4jl . On 
the other hand, the distance between hard spheres in contact must be equal 
to the diameter a = 1 of the spheres: 

11^-^-11 = 1 (2.6) 

Eq. ()2.6|) brings exactly N c constraints on the positions of the centers of the 
particles. Once again, there are Nd — d(d + l)/2 degrees of freedom for the 
particle positions. Therefore one must have Nd — d(d + l)/2 > iV c which 
implies z < 2d (note that this argument is contradicted in the case of the 
crystal whose coordination is larger. In the crystal, the constraints on the 
particle positions are redundant. Adding an infinitesimal poly-dispersity 
destroys this effect). Finally these two bounds lead to z = 2d. 

The second remark is that an isostatic state is marginally rigid: if one 
contact is removed, one soft mode appears. This has the following conse- 
quence: an isostatic state is very sensitive to boundary conditions. Consider 
a subsystem of size L in a large isostatic system. Let N ext ~ L d ~ x be the 
number of contacts of this subsystem with external beads. The number of 
contacts inside the subsystem N int is on average (N int ) = Nd — N ext /2, where 
the factor | shows up because a contact is shared by two particles. This im- 
plies that if the N ext contacts were to be removed, the subsystem would not 
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be rigid. On average it will have N ext /2— d(d+l)/2 soft modes. Consider now 
the force field composed of the N ext contact forces applied by the external 
beads on the subsystem. It must be orthogonal to the N ext /2 — d(d+l)/2 soft 
modes. This implies that roughly half of the external contact forces are free 
degrees of freedom. If the contact forces are imposed on half of the bound- 
ary of the subsystem, the contact forces of the other half of the boundary 
are determined. This is very different from an elastic body where one can 
impose any stress on the whole boundary and compute the response of the 
system. In an isostatic system forces propagate from one side of the system 
to the other. In |57j . the authors discuss the nature of the soft modes of iso- 
static anisotropic system. Using Eq. (|2.5|) they infer a null-stress law among 
the components of the stress tensor. This leads to the hyperbolic equations 
proposed to describe force propagation in granular matter fTJ E] • 

2.3 Covalent glasses 

Phillips [55^ used the soft modes, or "floppy modes", to study the structure 
of covalent glasses. The counting of degrees of freedom is slightly different 
from a network of springs where only the stretching of the contacts matters. 
Covalent interactions also display multi-body forces: the energy of the system 
depends on the angles formed by the different covalent bonds of an atom. 
These extra-terms in the energy, the "bond bending energy", bring extra- 
constraints on the soft modes. How many total constraints are there per atom 
of valence vl As we discussed with Eq. fl2.3j) . the stretching of the v bonds 
leads tov/2 constrains on the soft modes (as there is one constraint per bond 
and each bond is shared by two atoms). Furthermore, the covalent bonds 
form v(v — l)/2 angles, each of them corresponds to a term in the energy 
expansion. Therefore the total number of constraints is v/2 + v(v — l)/2 = 
v 2 /2. 

Consider, following 55,, chalcogenide alloys such as Ge,j;Sei_ x ., where the 
relative concentration x of Ge can vary from to 1. Ge has a valence of 4 
whereas Se has a valence of 2. When x = the glass as a polymeric structure. 
When x increases the connectivity of the covalent network increases too, until 
it becomes rigid. This takes place when x ■ 4 2 + (1 — x) ■ 2 2 = 2d = 6, that is 
when x c ~ 0.16. 

However, although the covalent network is floppy for smaller x, the glass 
is still rigid at lower concentration. In particular the shear modulus does 
not vanish below x c jHSj- The reason is that this counting argument neglects 
the weaker interactions induced by the lone pair electrons, that lead to short 
range repulsion and long range attraction (Van der Waals interaction). Nev- 
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ertheless several experimental results show that the transition at x — 0.16 
affects certain properties of the system. In Ge-As-Se glasses the fragility, that 
quantifies how the dependence of the viscosity with temperature is different 
from an Arrhenius law, is maximal at the transition [HI]. Furthermore, it 
was argued that the composition of the best glass former should be x c , as 
it appears experimentally, see Fig (j2.2|) . A qualitative argument is as follows: 
when x increases toward x c , the covalent network becomes more and more 
intricate and the viscosity of the system increases. It takes therefore more 
and more time to nucleate the crystal. On the other hand, if x increases 
above x c , the covalent network is over-constrained, the bonds are frustrated 
and have to store energy. The configuration of the crystal, where the bonds 
can organize to avoid this frustration, becomes more and more favorable with 
x in comparison with the amorphous state, and is thus easier to nucleate. 

This description is "mean field" as it does not consider the possible spa- 
tial fluctuations of the coordination. Such fluctuations may lead to a system 
where rigid, high-coordinated regions coexist with floppy, weakly-coordinated 
regions. To study this possibility models were proposed such as rigidity per- 
colation, see e.g. (HI] and reference therein. In its simplest form, this model 
considers springs randomly deposited on a lattice. When the concentration 
of springs increases, there is a transition when a rigid cluster percolates in 
the system almost fully floppy. Such a cluster is a non-trivial fractal object 
and contains over-constrained regions. In this form, this model is at infinite 
temperature, as there are no correlation among the contacts deposited. It 
is interesting to note the difference with the transition of jamming that we 
study in what follows, which is also a transition of rigidity. At the jam- 
ming threshold the correlations among particles are obviously important, 
as the temperature is not infinite. The system is exactly isostatic, as we 
shall discuss, and does not contain over-constrained regions. It is a normal 
d-dimensional object with only very few holes of size of order unity 4 , the 
"rattlers" , which are isolated particles without contact [SD] • 

4 The absence of large voids, or rattlers, has geometrical origins. Because particles are 
repulsive, the rigid region surrounding a void must be convex. A large void necessitates the 
creation of a vault. A flat vault is impossible as forces cannot be balanced on any particle 
of the surface. A weakly curved vault imposes drastic constraints on the distribution of 
angles of contact between these particles. This suggests that the probability of making a 
large void decays very fast, presumably at least exponentially in the surface of the void. 
If friction is present such voids might form more easily. 
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Fig. 2.2: The quenching rate (or difficulty of glass formation) was plotted 
as a function of x in Ge x Sei_ x alloys. Both lines are sketched 
in to guide the reader's eye. The dashed line corresponds to the 
network effect discussed in the text. Phillips |HHj . 
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2.4 The rigidity of "soft solids" 

Alexander noticed [SSI ^he following contradiction: there are solids less con- 
nected than required by the Maxwell criterion, and that are nevertheless 
rigid. For example, one may describe a gel as an assembly of reticulated 
point linked by springs that are the polymers. In general the coordination 
of such system is less than 6. Why then a gel has a finite shear modulus? 
The starting point of Alexander's answer is the simple following remark: the 
stress influences the frequency of the vibrations. At a macroscopic level the 
presence of a negative stress (that is, when the system is stretched) is re- 
sponsible for the fast vibrations of strings or drumheads. When a system 
is compressed, the stress is positive and lowers the frequency of the modes 
that can even become unstable, as for the buckling of a thin rod. To discuss 
the role of stress at a microscopic level, consider particles interacting with a 
potential V(r). The expansion of the energy leads to: 

SE = £ V'{r%)d rij + \v{^)dn% + 0{dr%) (2.7) 

ij 

where the sum is over all pairs of particles, rfj is the equilibrium distance 
between particles % and j. In order to get an expansion in the displacement 
field 8E4 we use: 

dm = {8Rj - SRj.fiij + [^gj_LpLJ! + (8R 3 ) (2.8) 

Zr ij 

where (SRj — 5Ri) ± indicates the projection of SRj — SRi on the plane or- 
thogonal to Hij. When used in Eq. (j2.7|) . the linear term in the displacement 
field disappears (the system is at equilibrium) and we obtain: 

5E = {£ V'(rtJ) m \i Rl)±]2 } + \v"{r%)[{6R 3 - 5R).n^ + 0(5R 3 ) 

(2-9) 

The difference from Eq. (j2.1|) is the term inside curly brackets. Such term 
is called the "initial stress" term in |HH! as it is directly proportional to the 
forces V'{r1j). We shall also refer to it as the applied stress term. If the 
system has a negative pressure p this term increases the frequency of the 
modes. As a consequence if p < all the soft modes of a weakly-connected 
system gain a finite positive energy: the system is rigid. This occurs in gels: 
the osmotic pressure of the solvent is larger than the external pressure. This 
imposes that the network of reticulated polymers carry a negative pressure 
to compensate this difference: the polymers are stretched. This rigidities the 
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system. As a consequence, the shear modulus of gels is directly related to 
the osmotic pressure of the solvent. 



3. Vibrations of isostatic systems 



3.1 Isostaticity 

When a system of repulsive spheres jams at zero temperature, the system 
is isostatic |HZ1 EH EE] (when the rattlers, or particles without contacts, are 
removed). As we said above, on the one hand it must be rigid and satisfy 
the bound of (|2.4|) . One the other hand, it cannot be more connected: this 
would imply that the contacts are frustrated as they cannot satisfy Eq. (j2.6|) . 
That is, there would not be enough displacements degrees of freedom to 
allow particles in contact to touch each other without interpenetrating, as 
we discussed for hard spheres in the last Chapter. Thus the energy of the 
system, and the pressure, would not vanish at the transition. It must be 
the case, since an infinitesimal pressure can jam an athermal gas of elastic 
particles. 

In terms of energy expansion, since the pressure vanishes at the transition, 
the initial stress in bracket in Eq. (j2.9j) vanishes. For concreteness in the 
following Chapters we consider the harmonic potential, corresponding to a = 
2 in Eq. (|l.l|) . The expansion of the energy is then given by Eq. (|2.1|) . In 
Chapter [7| we generalize our findings to other soft sphere potentials and 
other types of interactions. 

An isostatic system is marginally stable: if q contacts are cut, a space 
of soft modes of dimension q appears. For our argument below we need 
to discuss the extended character of these modes. In general when only 
one contact (ij) is cut in an isostatic system, the corresponding soft mode 
is not localized near (ij). This comes from the non-locality of the isostatic 
condition that gives rise to the soft modes; and was confirmed in the isostatic 
simulations of Ref |57j . which observed that the amplitudes of the soft modes 
spread out over a nonzero fraction of the particles. This shall be proved by 
the calculation of Chapter [7| that shows that in an isostatic system, the 
response to a local strain does not decay with the distance from the source. 
When many contacts are severed, the extended character of the soft modes 
that appear depends on the geometry of the region being cut. If this region 
is compact many of the soft modes are localized. For example cutting all the 
contacts inside a sphere totally disconnects each inner particle. Most of the 
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Fig. 3.1: One soft mode in two dimensions for iV « 1000 particles. The rela- 
tive displacement of the soft mode is represented by a line segment 
extending from the dot. The mode was created from a previously 
prepared isostatic configuration, periodic in both directions, fol- 
lowing [HD] . 20 contacts along the vertical edges were then removed 
and the soften modes determined. The mode pictured here is an 
arbitrary linear combination of these modes. 

soft modes are then the individual translations of these particles and are not 
extended throughout the system. 

In what follows we will be particularly interested in the case where the 
region of the cut is a hyper-plane as illustrated in Fig. ()3.2|) . In this situation 
occasionally particles in the vicinity of the hyper-plane can be left with less 
than d contacts, so that trivial localized soft modes can also appear. However 
this represents only a finite fraction of the soft modes. We expect that there 
is a non-vanishing fraction q' of the total soft modes that are not localized 
near the hyper-plane. Rather, as when a single contact is cut, these modes 
should extend over the whole system, like the mode shown in Fig. (|3.1|) . We 
shall define extended modes more precisely in the next section. 
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3.2 Variational procedure 

We aim to show first that the density of states of an isostatic system does not 
vanish at zero frequency. Since D(u) is the total number of modes per unit 
volume per unit frequency range, we have to show that there are at least of 
the order of ujL d normal modes with frequencies smaller than u for any small 
uj. As we justify later, if proven in a system of size L for uj ~ ujl ~ 1/A this 
property can be extended to a larger range of uj independent of L. Therefore 
it is sufficient to show that there are of the order of L d ~ l normal modes 
with frequency of the order of 1/L, instead of the order of one such mode 
in a continuous solid: the whole translation of the system. To do so we use 
a variational argument: M. is a positive symmetric matrix. Therefore if a 
normalized mode has an energy SE, we know that the lowest eigenmode has 
a frequency ujq = y/~E~o < \f5E. Such argument can be extended to a set 
of modes 5 : if there are m orthonormal trial modes with energy 5E < uj 2 , 
then there are at least m/2 normal modes with frequency smaller than \f2uj t . 
Therefore we are led to find of the order of L d ~ x trial orthonormal modes 
with energy of order 1/L 2 . 

3.3 Trial modes 

For concreteness we consider the three-dimensional cubic iV-particle system 
S of Ref [SHj with periodic boundary conditions at the jamming threshold. 
We label the axes of the cube by x, y, z. S is isostatic, so that the removal 
of n contacts allows exactly n displacement modes with no restoring force. 
Consider for example the system S' built from S by removing the q ~ L 2 
contacts crossing an arbitrary plane orthogonal to (ox); by convention at 
x = 0, see Fig. (|H.2|) . S', which has a free boundary condition instead of 
periodic ones along (ox), contains a space of soft modes of dimension q 6 , 
instead of one such mode — the translation of the whole system — in a normal 
solid. As stated above, we suppose that a subspace of dimension q' ~ L 2 of 
these soft modes contains only extended modes. We define the extension 
of a mode relative to the cut hyper-plane in terms of the amplitudes of 
the mode at distance x from this hyper-plane. Specifically the extension 

5 If m a is the a'th lowest eigenvalue of M. and if e a is an orthonormal basis such that 
{e a \M.\e a ) = n a then the variational bound of A. Horn [Am. J. Math 76 620 (1954)] shows 
that J2i m a < J2i n a- Since qn q > J2i n a> ancl since J2l m a > H\l2 m a > (<l/2) m q/2, 
we have qn q > (q/2)m q / 2 as claimed. 

6 The balance of force can be satisfied in 5' by imposing external forces on the free 
boundary. This adds a linear term in the energy expansion that does not affect the 
normal modes. 
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e of a normalized mode |5R) is defined by J2i sin 2 (^)(2|<5R) 2 = e, where 
the notation (i\5TV) indicates the displacement of the particle i of the mode 
considered. For example, a uniform mode with (i\5TV) constant for all sites 
has e = | independent of L. On the other hand, if (i|<5R) = except for a site 
i adjacent to the cut hyper-plane, the Xi/L ~ L^ 1 and e ~ L~ 2 . We define 
the subspace of extended modes by setting a fixed threshold of extension e 
of order 1 and thus including only soft modes (3 for which ep > e$. As we 
argued at the beginning of the Chapter, we expect that a fraction of the total 
soft modes are extended. Thus if q' is the dimension of the extended modes 
vector space, we shall suppose that q'/q remains finite as L — > oo; i.e. a fixed 
fraction of the soft modes remain extended as the system becomes large. The 
appendix at the end of this Chapter presents our numerical evidence for this 
behavior. 

We use them to build q' orthonormal trial modes of frequency of the order 
1/L in the initial system S. Let us denote |5R/j) a normalized basis of the 
vector space of such extended modes, 1 < /3 < q'. These modes are not soft 
in the jammed system S because they deform the previous q contacts located 
at x = 0, and therefore cost energy. Nevertheless a set of trial modes, |5RJ), 
can still be formed by altering the soft modes so that they do not have an 
appreciable amplitude at the boundary where the contacts were severed. We 
seek to alter the soft mode to minimize the distortion at the severed contacts 
while minimizing the distortion elsewhere. Accordingly, for each soft mode 
j3 we define the corresponding trial-mode displacement (i|<5R*) to be: 

(i\5R;) = C p sm(^)(i\5Rp) (3.1) 

where the normalization constant Cp depends on the spatial distribution of 

the mode (3. If for example (i|5R) = except for a site i adjacent to the cut 

plane, Cp grows without bound as L — >• oo. In the case of extended modes 

Cp 2 = s i n2 (^f E )0I^R/3) 2 — e /3 > e o, and therefore Cp is bounded above 
_i 

by e 2 . The sine factor suppresses the problematic gaps and overlaps at the 
q contacts near x = and x = L. Formally, the modulation by a sine is a 
linear mapping. This mapping is invertible if it is restricted to the extended 
soft modes. Consequently the basis |5R/?) can always be chosen such that 
the |<5Rjg) are orthogonal. Furthermore one readily verifies that the |5R^)'s 
energies are small, because the sine modulation generates an energy of order 
1/L 2 as expected. Indeed we have from Eq. (j2.1|) : 

SE = C 2 J2[( S m(^)(i\5Rp) - sm(^)(j\5Rp}) ■ r%] 2 (3.2) 

Using Eq. ()2.3|) . and expanding the sine, one obtains: 
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Fig. 3.2: Illustration of the boundary contact removal process described in 
the text. Eighteen particles are confined in a square box of side 
L periodically continued horizontally and vertically. An isostatic 
packing requires 33 contacts in this two-dimensional system. An 
arbitrarily drawn vertical line divides the system. A contact is 
removed wherever the line separates the contact from the center of 
a particle. Twenty-eight small triangles mark the intact contacts; 
removed contacts are shown by the five white circles. 
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5E « ^£cos 2 (^)^(n^ ■ <Q 2 «j\6Rp) ■ n^f (3.3) 

(ij) 

<e \n/Lfj:(j\^) 2 (3-4) 

where e x is the unit vector along (ox). The sum on the contacts can be 
written as a sum on all the particles since only one index is present in each 
term. Using the normalization of the mode (3 and the fact that the coordi- 
nation number of a sphere is bounded by a constant z max (z max = 12 for 3 
dimensional spheres 7 ), one obtains: 

SE < e Q l {Ti/L) 2 z max = u\ (3.5) 

One may ask if the present variational argument can be improved, for exam- 
ple by considering geometries of broken contacts different from the surface 
we considered up to now. When contacts are cut to create a vector space of 
extended soft modes, the soft modes must be modulated with a function that 
vanishes where the contacts are broken in order to obtain trial modes of low 
energy. One the one hand, cutting many contacts increases the number of 
trial modes. On the other hand, if too many contacts are broken, the mod- 
ulating function must have many "nodes" where it vanishes. Consequently 
this function displays larger gradients and the energies of the trial modes 
increase. Cutting a surface (or many surfaces, as we shall discuss below) ap- 
pears to be the best compromise between these two opposing effects. Thus 
our argument gives a natural limit to the number of low-frequency states to 
be expected. 

Finally we have found of the order of L 2 trial orthonormal modes of 
frequency bounded by tu^ ~ and we can apply the variational argument 
mentioned above: the density of statesis bounded below by a constant below 
frequencies of the order ujl- In what follows, the trial modes introduced 
in Eq. (j3.1)) . which are the soft modes modulated by a sine, shall be called 
"anomalous modes" . 

7 In a polydisperse system z max could a priori be larger. Nevertheless Eq. I|3.3|I is a sum 
on every contact where the displacement of only one of the two particles appears in each 
term of the sum. The corresponding particle can be chosen arbitrarily. Chose the smallest 
particle of each contact. Thus when this sum on every contact is written as a sum on 
every particles to obtain Ea. (|3.5J) . the constant z max still corresponds to the monodisperse 
case, as a particle cannot have more contacts with particles larger than itself. 
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3.4 Argument extension to a wider 
frequency range 

We may extend this argument to show that the bound on the initial density 
of states extends to a plateau encompassing a nonzero fraction of the modes 
in the system. If the cubic simulation box were now divided into m 3 sub- 
cubes of size L/m, each sub-cube must have a density of states equal to the 
same D{uj) as was derived above, but extending to frequencies of order mui- 
These subsystem modes must be present in the full system as well, therefore 
the bound on D(u) extends to [0, itiujl]. We thus prove that the same bound 
on the average density holds down to sizes of the order of a few particles, 
corresponding to frequencies independent of L. Finally D{ui) does not vanish 
when uj — > 0, as indicates the presence of the observed plateau in the density 
of states. We note that in d dimensions this argument may be repeated to 
yield a total number of modes, below a frequency ujl ~ thus 

yielding a limiting nonzero density of states in any dimension. 

We note that the trial modes of energy 5E ~ Z -1 that we introduced 
by cutting in subsystems of size / are, by construction, localized to distance 
scale /. Nevertheless we expect these trial modes to hybridize with the trial 
modes of other subsystems, and the corresponding normal modes not to be 
localized on such length scale. 

3.5 Appendix: Spatial distribution of the 
soft modes 

In our argument we have assumed that when q ~ L^ 1 contacts are cut 
along a hyper-plane in an isostatic system, there is a vector space of di- 
mension q' = aq which contains only extended modes, such that a does not 
vanish when L — > oo. A normalized mode \8H) was said to be extended if 
Y^i sin 2 (i\5TV) 2 > eo, where eo is a constant, and does not depend on L. 
Here we show how to choose eo > so that there is a non-vanishing fraction 
of extended soft modes. We build the vector space of extended soft modes 
and furnish a bound on its dimension. 

Let us consider the linear mapping Q which assigns to a displacement 
field |5R) the displacement field (i\Q5TL) = sin 2 (xj7r/L)5i?j. For any soft 
mode |5R#) one can consider the positive number ap = (^R^C^R^) = 
2 -* 2 

X^iSin (xiTi/ L)5R/3 i . We build the vector space of extended modes by recur- 
sion: at each step we compute the dp for the normalized soft modes, and we 
eliminate the soft mode with the minimum dp. We then repeat this procedure 
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Fig. 3.3: The overlap function f(x) as defined in Eq. ()3.6j) for different sys- 
tem sizes in three dimensions. Soft modes were created from iso- 
static configurations as described below Fig. 13.11 

in the vector space orthogonal to the soft modes eliminated. We stop the 
procedure when cjg > e for all the soft modes (3 left. Then all the modes left 
are extended according to our definition. We just have to show that one can 
choose e > such that when this procedure stops, there are q' modes left, 
with q' > aq and a > 0. In order to show that, we introduce the following 
overlap function: 



The sum is taken on an orthonormal basis of soft modes (3 and on all 
the particles whose position has a coordinate x,- t G [x,x + dx\. f(x) is the 
trace of a projection, and is therefore independent of the orthonormal basis 
considered. f(x) describes the spatial distribution of the amplitude of the 
soft modes. The are normalized and therefore: 



We have examined soft modes made from configurations at the jamming 
transition found numerically in |5()j . The overlap function f(x) was then 
computed for different system sizes L. These are shown in Fig. 13.31 It 
appears from Fig. ()3.3|) that i) when f(x) is rescaled with the system size 
it collapses to a unique curve, and ii) this curve is bounded from below 
by a constant c (c ~ 0.6). Consequently one can bound the trace of Q: 
trQ = Jq qf(x) sin 2 (x) > qc/2. On the other hand one has trQ = J2p=i a i3, 
where the sum is made on the orthonormal basis we just built. Introducing 



f(x)dx^q- 1 J2 E Md 2 



(3.6) 



j3=l,...,q Xi£[x,x+dx] 




(3.7) 
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the numbers a and e such that there are q' = aq extended modes, and using 
that the ap < 1, one can bound this sum and obtain trQ < aq + (1 — a)qeo- 
Choosing for example eo = c/2 > 0, one finds a > c/(2 — c), which is a 
constant independent of L as claimed. 



4. Evolution of the modes with the 
coordination 



4.1 An "isostatic" length scale 

When the system is compressed and moves away from the jamming transi- 
tion, the simulations showed that the extra-coordination number Sz = z — z c 
increases. In the simulation, the compression also creates forces on all con- 
tacts. In this Chapter we ignore these forces, and instead only consider the 
contact network created by compression, but in the absence of applied pres- 
sure. Any tension or compression in the contacts is removed. The effect on 
the energy is to remove the first bracketed term from Eq. (j2.9j) above, and 
the expansion of the energy is still given by Eq. (|2.1|) . We note that removing 
these forces, which add to zero on each particle, does not disturb the equi- 
librium of the particles or create displacements. In this section we ignore the 
question of how 5z depends on the degree of compression. We return to this 
question in the next section. Compression causes AN C = N5z/2 ~ L 3 5z ex- 
tra constraints to appear in Eq. (j2.1j) . Cutting the boundaries of the system, 
as we did above, relaxes q ~ L 2 constraints. For a large system, q < AN C 
and Eq. (j2.1|) is still over-constrained so that no soft modes appear in the sys- 
tem. However, as the systems become smaller, the excess AN C diminishes, 
and for L smaller than some I* ~ 5z~ l where AN C = q, the system is again 
under-constrained, as was already noticed in [HZj- This allows us to build 
low-frequency modes in subsystems smaller than I*. These modes appear 
above a cut-off frequency u* ~ Z* -1 ; they are the excess-modes that con- 
tribute to the plateau in D(uS) above to*. In other words, anomalous modes 
with characteristic length smaller than I* are little affected by the extra con- 
tacts, and the density of states is unperturbed above a frequency oj* ~ Sz. 
This scaling is checked numerically in Fig ]4.1l It is in very good agreement 
with our prediction up to Sz ~ 2. 

At frequency lower than uo* we expect the system to behave as a disor- 
dered, but not ill-connected, elastic medium, so that the vibrational modes 
are similar to the plane waves of a continuous elastic body. We refer to these 
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Fig. 4.1: uj* as defined in the insert of Fig.l vs. coordination, in the system 
with relaxed springs. The line as a slope one. 

modes as "acoustic modes". Thus we expect D(uj) at small uj to vary as 
uj 2 c~ 3 , where c(5z) is the sound speed at the given compression. This c may 
be inferred from the bulk and shear moduli measured in the simulations; that 
we shall derive in chapter One finds for the transverse velocity q ~ (5z)z, 
and for the longitudinal velocity q ~ 5z°. Thus at low frequency the D(oS) is 
dominated by the transverse plane waves and at uj = uj* the acoustic density 
of states is uj 2 c 3 ~ 8z 2 8z~ 3 / 2 ~ 5z^\ the acoustic density of states should be 
dramatically smaller than the plateau density of states. There is no smooth 
connection between the two regimes, thus we expect a sharp drop-off in D(uj) 
for to < uj* . Such drop-off is indeed observed, as seen in Fig (j4.2j) . In fact, 
because of the finite size of the simulation, no acoustic modes are apparent 
at uj < uj* near the transition. 

Thus the behavior of such system near the jamming threshold depends on 
the frequency uj at which it is considered. For uj > uj* the system behaves as 
an isostatic state: the density of states is dominated by anomalous modes. 
For uj < uj* we expect it to behave as a continuous elastic medium with 
acoustic modes. Since the transverse and the longitudinal velocities do not 
scale in the same way, the wavelengths of the longitudinal and transverse 
plane waves at uj* are two distinct length scales k and It which follow k ~ 
quj*" 1 and l t ~ c t uj*~ l . At shorter wave lengths we expect the acoustic 
modes to be strongly perturbed. Note that since q ~ 5z°, one has li ~ I*. 
Interestingly, l t ~ 5z~z is the smallest system size at which plane waves can 
be observed: for smaller systems, the lowest frequency mode is not a plane 
wave, but an anomalous mode, k was observed numerically in 
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Fig. 4.2: Log-linear plot of the density of states for N=1024 for three values 
of <p — <f) c in the soft spheres system (dotted line) and the system 
where the applied stress term has been removed (solid line). 

4.2 Role of spatial fluctuations of z 

Our argument ignores the spatial fluctuations of 5z. If these fluctuations 
were spatially uncorrelated they would be Gaussian upon coarse-graining: 
then the extra number of contacts AN C in a subregion of size L would have 
fluctuations of order L d l 2 . The scaling of the contact number that appears 
in our description is AiV c ~ L d ~ l and is therefore larger than these Gaussian 
fluctuations for d > 2. In other terms at the length scale I* where soft 
modes appear, the fluctuations of the number of contacts inside the bulk 
are negligible in comparison with the number of contacts at the surface. 
Therefore the extended soft modes that are described here are not sensitive 
to fluctuations of coordination in three dimensions near the transition. In 
|fi7j we argued that in two dimensions there are spatial anti-correlations in 
z, and that the fluctuations do not affect the extended soft modes in two 
dimensions either. 

Note that these arguments do not preclude the existence of low-frequency 
localized modes that may appear in regions of small size i < 1*, and that 
could be induced by very weak local coordination, or specific configuration. 
The presence of such modes would increase the density of states at low- 
frequency. There is no evidence for their presence in the simulations of [50] . 
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Fig. 4.3: Rigid unit modes model applied to silica. Trachenko et al. 68J. 

4.3 Application to tetrahedral networks 

The nature of slow vibrations -and the possible presence of two-level systems- 
has been much studied in silica, one of the main glass-forming system. In this 
glass (or more generally aluminosilicates) the forces within the tetrahedra 
S2O4 are much stronger than the forces that act between them |69j : it is 
easier to rotate two linked tetrahedra than to distort one tetrahedron 8 . This 
suggests to model such glass as an assembly of linked tetrahedra loosely 
connected at corners: this is the "rigid unit modes" model [ZOj- In this 
model the tetrahedra are characterized by a unique parameter, a stiffness k 
9 . Recently this model was used to study the vibrations of silica [71]. The 
authors first generate realistic configurations of Si0 2 at different pressures 
using molecular dynamics simulations. At low pressure, they obtain a perfect 
tetrahedral network. When the pressure becomes large, the coordination 
of the system increases with the formation of 5-fold defects. Once these 
microscopic configurations are obtained, the rigid unit model is used and the 
system is modeled as an assembly of rigid tetrahedra, see Fig. (j4.3J) . Then, 
the density of states of such network is computed. The results are shown 

8 For example the bending energy of Si-o-Si is roughly 10 times smaller than the stretch- 
ing of the contact Si-o [73] . 

9 In fact the rigidity of a tetrahedron induced by the covalent bonds should be char- 
acterized by 3 parameters corresponding to different deformations of the tetrahedron. If 
these parameters are of similar magnitude, as one expects for example for silica, this does 
not change qualitatively the results discussed here. 
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Fig. 4.4: Density of rigid unit modes for silicate at different pressure. Tra- 
chenko et al.|7Tj. 



in Fig. (|4.4jl . One can note the obvious similarity with the density of states 
near jamming of Fig. (jl.3|) . We argue that the cause is identical, and that 
the excess-modes correspond to the anomalous modes made from the soft 
modes, rather than to one-dimensional modes as proposed in [71]. Indeed, 
a tetrahedral network is isostatic, see e.g. [HE]- The counting of degrees of 
freedom can be made as follows: on the one hand each tetrahedron has 6 
degrees of freedom (3 rotations and 3 translations). On the other hand, the 
4 corners of a tetrahedron bring each 3 constraints shared by 2 tetrahedra, 
leading to 6 constraints per tetrahedron. Thus the system is isostatic. When 
the pressure increases the coordination increases too, leading to the erosion of 
the plateau in the density of states discussed earlier in this Chapter. Finally, 
these predictions of the density of states fail to describe silica glass vibrations 
at low frequencies, where one cannot neglect the weaker interactions anymore 
nor the role of the initial stress that we discuss in the next Chapter. In 
Chapter |H] we evaluate the effect of the weaker interactions. This allows us 
to propose an explanation for the nature of the boson peak in such glasses. 



5. Effect of the initial stress on 
vibrations 



In this section we describe how the above simple description of D(u) is af- 
fected by the presence of applied stress. In general when a system of particles 
at equilibrium is formed, there are forces between interacting particles. For 
harmonic soft spheres it leads to a non- vanishing first term in Eq. (j2.9j) that 
becomes: 



where we used « 1. This term in bracket is (a) negative for repulsive 
particles (b) proportional to the transverse relative displacement between 
particle in contact (c) scales as the pressure p, and is therefore vanishing at 
the jamming transition. The full dynamical matrix T> can be written: 



where M 1 is written in tensorial notation in footnote 10 . The spectrum of 
T> has a priori no simple relation with the spectrum of M.. Because M! is 
much smaller than M. near the transition, one can successfully use pertur- 
bation theory for the bulk part of the normal modes of M.. Nevertheless 
perturbation theory fails at very low frequency, which is of most interest. 
In this region the spectrum of M. contains the plane waves and the anoma- 
lous modes. In what follows we estimate the change of frequency induced 
by the applied stress on these modes. We show that the relative correction 
to the plane wave frequencies is very small, whereas the frequency of the 
anomalous modes can be appreciably changed. Finally we show that these 
considerations lead to a correction of the Maxwell criterion of rigidity. 

10 in three dimension we have Ai'ij — ~ ^r^' 3 [S(ij){fhij®mij+kij($ikij)+6i t j J2<i>( f ™ij® 
rhij + kij <8> kij), where (n,j ,rhij ,kij) is an orthonormal basis. 




V = M+M' 
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5.1 Applied stress and plane waves 

Consider a plane wave of wave vector k. Since the directions random, 
both the relative longitudinal and the transverse displacements are of the 
same order: [(8B4 — SRj) 1 -] 2 ~ [{8B4 — 5Rj).nij} 2 ~ (5Ri) 2 k 2 . Consequently 
the relative correction AE/E induced by the applied stress term is very 
small: 

E ~ T, {ij) [{8R j -8R i )-n ij \ 2 

since rf? — 1 is proportional to the pressure p, while the others factors remain 
constant as p — > 0, ^ ~ p, and is thus arbitrarily small near the jamming 
threshold 11 . Note that when the pressure is high, this effect is non negligi- 
ble. In particular elastic instabilities can occur, and can be responsible for 
conformational changes, see [73] for such examples in silica crystals. 



5.2 Applied stress and anomalous modes 

For anomalous modes the situation is very different: we expect the transverse 
relative displacements to be much larger than the longitudinal ones. Indeed 
soft modes were built by imposing zero longitudinal terms, but there were no 
constraints on the transverse terms. These are the degrees of freedom that 
generate the large number of soft modes. The simplest assumption is that 
the relative transverse displacements are of the order of the displacements 

themselves, that is J2(ij)[(8Rj — SRi)- 1 ] 2 ~ J2i$Ri — 1 f° r the anomalous 
modes that appear above u*. This estimate can be checked numerically for an 
isostatic system where the sum of the relative transverse terms is computed 
for all uj. The sum converges to a constant when u — > as assumed, see 

Fig. dHm. 

We can estimate the scaling of the correction in the energy AE induced 
by the stress term on the anomalous modes: 

11 In disordered systems the acoustic modes are not exact plane waves, see e.g. the recent 
simulations in Lennard- Jones systems |B10IH1III|- As we discuss below, for transverse plane 
waves the energy is reduced by a factor Sz. Therefore the relative correction of energy 
induced by the applied pressure is of the order of rather than p. We have, as shown, 
j- <C 1, so that we still expect the correction to be small near the jamming threshold. In 
principle non-afhne displacements could have other interesting effects, such as an increase 
of the transverse terms amplitude. If so, the effect of applied stress on acoustic modes 
would be enhanced. 
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Fig. 5.1: Sum of the transverse terms (red curve) e = | J2(ij)[(8Rj ~ bRi) 1 '} 2 

and longitudinal terms (black curve) e = \ Y^(ij)[{o~Rj — 8Ri) • ^j] 2 
for each mode of frequency uo at the jamming threshold in three 
dimensions. The longitudinal term is equal to the energy of the 
modes and vanish quadratically at frequency. The transverse 
term converges toward a constant different from 0. 



AE --£(!- rj)[(^ - 5R^f ~ -p (5.4) 
to'} 

which is an absolute correction, which can be non-negligible in comparison 
with the energy E. 

5.3 Onset of appearance of the anomalous 
modes 

We can now estimate the lowest frequency of the anomalous modes. The 
modes that appear at u* in the relaxed springs system have an energy lowered 
by an amount of order — p in the original system. Applying the variational 
theorem of the last section to the collection of slow modes near uo* indicates 
that there must be slow normal modes with a lower energy. That is, the 
frequency u>am at which anomalous modes appear verifies: 

uam < [(uj*) 2 - A 2 p]^ = \A x bz 2 - A 2 p]^ (5.5) 

where A\ and A 2 are two positive constants. This indicates that the im- 
portant parameters of the low frequency excitations are coordination and 
stress. 
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Fig. 5.2: Phase diagram of rigidity in terms of coordination and pressure. 

When p > 0, the line separating the stable and unstable regions is 
defined by Eq. (j5.6|) . When the pressure is negative, any connected 
system can be rigid, as it is the case for gels. 

5.4 Extended Maxwell criterion 

From this estimation we can readily obtain a relation among coordination 
and pressure that guarantees the stability of a system. There should be no 
negative frequencies in a stable system, therefore ujam > 0. Thus in an 
harmonic system the right hand side of Eq. (j5.5|) must be positive: 

5z > C p^ = 5z min (5.6) 

where Cq is a constant. This inequality, which must hold for any spatial 
dimension, indicates how a system must be connected to counterbalance the 
destabilizing effect of the pressure. A phase diagram of rigidity is represented 
in Fig. (j5.2j) . When p = 0, the minimal coordination z c corresponds to the 
isostatic state: this is the Maxwell criterion. As we said earlier for spherical 
particles z c = 2d. As we discuss later z c = d + 1 for particles with friction. 
When p > 0, Eq. (j5.6|) delimits the region of rigid systems: for example 
granular matter or emulsions lie above this line. When p < 0, even systems 
far less connected than what the Maxwell criterion requires are rigid [35] , 
as we discussed in Chapter El These systems contain many soft modes as 
defined in Eq. (|2.Hjl . but there are all stabilized by the positive bracketed term 
of Eq. ()2.9|) . Note that a similar phase diagram, with the same singularity of 
5z but a different z c , was recently obtained by a mean field approach [71] . 

Finally the relation ()5.6|) is verified in the simulations of jSU] where 5z ~ 
p2. In the next Chapter we justify why the inequality ()5.6|) is in fact an 
equality in the system of [5U] and discuss what determines the microscopic 
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structure of this system. 



6. Microscopic structure and marginal 
stability 



In the previous section we studied how the low-frequency vibrational prop- 
erties were related to the microscopic structure. The applied pressure has 
two antagonist effects: on the one hand, it increases the coordination num- 
ber, which stabilizes the system and increase the frequency of appearance 
of the anomalous modes. On the other hand, the applied pressure appears 
in the expansion of the energy and lowers the frequency of the anomalous 
modes. In what follows we study the relative amplitude of these two effects, 
or equivalently, where such amorphous solids are located in the (z,p) plane 
of Fig. Q . 

Since these systems are out-of-equilibrium, their microscopic structures 
depend on the system history. As we shall see below, the preparation of 
the system of [20] leads to marginally stable systems at any 0. The two 
antagonists effects of the pressure compensate 12 , which leads to u* ^> ujami 
as it appears in Fig. (|4.2jl . and 5z ~ pa. In what follows we propose a 
simple argument to justify such a behavior. We discuss in particular (i) the 
dynamic that takes place when a liquid of repulsive spheres is hyperquenched 
and (ii) the decompression of a jammed solid at zero temperature. This will 
also enable us to discuss the surprising geometrical property of the random 
close packing evoked in the introduction: there is a divergence in the pair 
correlation function g(r) ~ 1 j \Jr — 1 at close contact. We propose that this 
divergence is a vestige of the marginal stability that occurs at higher packing 
fraction. 

6.1 Infinite quench 

The simulations of [SU] show 5z ~ p^, thus saturating the bound of Eq. ()5.6|) . 
so that there are excess modes extending to frequencies much less than to*. 

12 Assuming an exact compensation of these two terms lead to loam — in an infinite 
size system. In Fig. (j4.2|l losm is slightly different from zero, as one would expect for a 
finite size system. 
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We start by furnishing an example of dynamics that lead to such a situation. 
Consider an initial condition where forces are roughly balanced on every 
particle, but such that the inequality (|5.fij) is not satisfied. Consequently, this 
system is not stable: infinitesimal fluctuations make the system relax with 
the collapse of unstable modes. Such dynamics was described by Alexander in 
PU] as structural buckling events: they are are induced by a positive stress as 
for the buckling of a rod, but take place in the bulk of an amorphous solid. 
These events a priori create both new contacts and decrease the pressure. 
When the bound of (|5.fi|) is reached, there are no more unstable modes. 
If the temperature is zero, the dynamics stops. Consequently one obtains 
a system where Eq. (|5.6|) is an equality, therefore (i) this system is weakly 
connected (ii) ujam = 0, there are anomalous modes much different from 
plane waves extending to zero frequency. A similar argument is present in 

In the simulations of Ref [50J the relaxation precedes as follows. The 
system is initially in equilibrium at a high temperature. Then it is hyper- 
quenched to zero temperature. At short time scales the dynamics that fol- 
lows is dominated by the relaxation of the stable, high frequency modes. The 
main effect is to restore approximately force balance on every particle. At 
this point, if inequality (|5.fij) is satisfied, the dynamics stops. If it is not sat- 
isfied, we are in the situation of buckling described above. The pressure and 
co-ordination number continue to change until the last unstable mode has 
been stabilized. At this point the bound of Eq. (20) is marginally satisfied, 
and there is no driving force for further relaxation. 

It is interesting to discuss further which situations lead to marginal stabil- 
ity. We expect that the situation of marginal stability that follows an infinite 
cooling rate takes place for a domain of the parameters of initial conditions 
(4>, T) , located at high temperature and low density. This domain could end 
at a finite <fi even when temperature is infinite, as it does in some related 
systems. This was shown in simulations and theoretical work on Euclidian 
random matrix [42J were most of the unstable modes vanish beyond a finite 
even when T = oo. When the cooling rate is finite, as in experiments, 
we expect that the relaxation does not stop when all the modes are stable, 
but that there are activated events that lead to further collapses of anoma- 
lous modes. These events a priori increase the connectivity and decrease the 
pressure beyond the bound of Eq. (j5.6J) . leading to uam > 0. In agreement 
with this idea, hyperquenched mineral glasses show a much stronger excess 
of modes [7S] in comparison with normally cooled glasses, and annealed poly- 
meric glasses expectedly show a smaller excess of modes |77j . 
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6.2 Decompression 

The data of (HOI were obtained by gradually decreasing the pressure from 
the above initial state of zero temperature and nonzero pressure. When 
pressure is lowered, it is observed that the system remains marginally stable. 
Here we propose a qualitative interpretation of these findings in the case 
of harmonic particles 13 . In [201, the decompression is obtained by discrete 
steps. At each step, the radii of each particle is reduced by a small amount e, 
while the center of the particles are kept fixed. This corresponds to an affine 
decompression. This new configuration is not at equilibrium in general. Then 
the particles are let to relax. The affine decompression has two effects: on the 
one hand it causes some contacts to open, on the other hand it reduces the 
pressure. The opening of these contacts tends to destabilize normal modes 
and reduce their frequencies, while the reduction in pressure tends to stabilize 
them. As we argue below, the destabilizing effect dominates. Thus, the 
affine decompression leads the system into the unstable region. Therefore we 
expect that when the particles relax, one recovers the dynamics that follows 
an infinite quench: modes buckle. As for the infinite quench of temperature 
discussed above, buckling should occur as long as the relaxation of the stable 
normal modes, which is faster than the collapse of unstable modes, does not 
bring back the system into the stable region. If so, the buckling increases 
the contact number and decreases the pressure until marginal stability is 
achieved, so that the inequality of Eq. ()5.fi|) is marginally satisfied as the 
pressure decreases, as observed in the simulations of [HQ] , 

In what follows we justify the claim that the destabilizing effect of the 
opening of the contacts dominate the effect of the pressure reduction. When 
the particles radii decrease by an amount e, a certain fraction e ~ g(l)e 
of contacts opens, where g{r) denotes the radial distribution function. For 
harmonic particles we expect g(l) ~ {(f) — (pc)' 1 14 - Hence, using that p ~ 
{(f) — (f> c ) for harmonic particles jJU] -as we shall demonstrate in the next 
Chapter-, one obtains e ~ |. On the other hand, the affine decompression 
lowers the pressure by an amount 5p ~ 5(f) ~ e. Thus the system can only 

13 To extend this argument to other potentials, for example Hertzian contacts, one should 
assume that g(l) scales at least as (0 — <^ c ) _3 when the jamming threshold is approached. 
This would imply that the number of contacts lost when the system is decompressed is 
large enough to generate buckling. This point is related to the evolution of the force 
distribution of Hertzian particles near the isostatic point. It is a subtle issue, and we are 
not aware of any numerical results of this sort. 

14 This is related to well-known empirical facts of the force distribution: one has 
P{F)dF ~ g(r)dr. For harmonic particles dr ~ dF and therefore P{F) ~ g{r). When 
rescaled by ~ p, P{F) converges to a master curve with P(F/(F) = 0) ^ 15 1| . This 
implies that g(l) ~ p' 1 ~ {(f)- </> c ) _1 [7S|. 
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afford to lose a fraction / of contacts while remaining stable: according to 

Eq. (15.61) : / = d ^ 5z " lin > Sp ~ p~^5p ~ 4 < e. Therefore / < e as claimed. 

a P pi 

Hence if an affine reduction of packing fraction e is imposed, far too many 
contacts open and the system is unstable. 

To conclude, it was observed in the simulations of (SI] that if the steps 
e are small, the decompression that takes place in jSHl is reversible: cycles 
of decompression-compression bring the system back to its initial configura- 
tion. This empirical fact indicates the absence of discontinuous irreversible 
events. Thus, the buckling generated by the opening of few contacts when 
e is small enough does not lead to rearrangements of finite amplitude much 
larger than e. This indicates that the dynamic of modes collapse increases 
the coordination by re-closing most of the contacts that open during the 
affine decompression (whose particles are separated by distance of at most 
e), and not by forming new contacts. It is reasonable to think that if several 
cycles of compression/decompression are made, the system will end up to be 
reversible. Why it is already so at the first decompression is a subtle question 
that we do not try to justify here. 

6.3 g(r) at the random close packing 

The probability g{r) of having two particles separated by a vector of length r 
displays a square root divergence g(r) ~ (r — at the jamming transition. 
In the introduction we pointed out that this divergence corresponds to the 
singular increase of the excess coordination 5z with the packing fraction, 
as was noted in jJU]. Here justify further this correspondence. We show 
that, if the decompression is assumed to be adiabatic, the singularity in g(r) 
is a necessary consequence of the marginal stability that characterizes the 
decompression. We argue that the pair of particles almost touching at C , 
responsible for the divergence in g(r), are the vestiges of the contacts that 
were present at higher to stabilize the system. In order to show this, we 
first have to count the contacts that open for a given — <f> c . Then we shall 
estimate the distance between the corresponding pairs of particles at the 
jamming threshold. 

As we discussed in the last section, when the system is decompressed 
it remains marginally stable: the coordination follows Eq. (j5.6|) . and z = 
z c + 5z min . Thus the density of contacts n(<p) per unit that open for a 
given follows exactly in the large L limit: 



n(<f>) 





1 
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We now would like to evaluate the distance that separates such particles 
at the jamming threshold. Let be the random variable that corresponds to 
the spacing r — 1 at the jamming threshold between two neighboring particles 
whose contact opened at a given (p. We want to estimate the fluctuations 
of Hq. If the decompression was purely affine would be single-valued and 
given by K ff {4>): 

K,M - (6.2) 

As we discussed in the last section the displacements that follow a decom- 
pression are not affine. For our present argument we need to evaluate the 
variance of the distribution of h^. It is directly related to the variance of the 
non-affine displacements that appear while decompressing. We expect that 
such non-affine displacements simply lead to a variance of of order of its 
average h a ff(4>) 15 . Therefore the probability P^,(h) that two particles whose 
contact opened at a given are at a distance 1 + h at the threshold can be 
written as: 

P *W = irr^-rhr^ (6 - 3) 

<P - <Pc <P- <Pc 

where / is a normalized scaling function / f(x)dx = 1. Thus one can compute 
g(r) by summing over all the contributions of the contacts that opened at 
> 4> c , as we represent on Fig. fjfi.lj) : 

g(r) ~ / n(4>)P^h = r - l)d<j) - 




(0_0 c )3/2^ 



1 f l 

9( r ) ~ 7 7TT / u-if(u)du (6.5) 

(r — 1) 2 J 

15 More generally, if two particles (in contact or not) are at a distance r of order one at 
a given <p, we expect that the fluctuations of the distance that separate them at </> c due to 
non-affine displacements is of order of the affine increase of their distance Sr = r ^~^ c w 



d4 ~ ^af f (0) ■ This comes from the following observation: non-affine displacements are 
induced by the requirement of stability that creates correlations among particles motions. 
To evaluate such correlations, consider the typical situation discussed in the last section 
where particles in contact at a given <fr have to stay in contact until <j> c , instead of spreading 
apart if the displacement was affine. At <f> c the inter-particle distance is 1, instead of a value 
r < 1 + h a f f (4>) . Thus we evaluate the typical departure from a pure affine displacement 
to be of order h a ff{<fi)- 
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h 




<t>c <t>l <t>2 



Fig. 6.1: Representation of typical variations of the spacing h between 



particles whose contacts opened at given > C . The arrows are 
typical trajectories. Those starting from 0! are thicker than those 
starting from 02, indicating that the density of contacts that open 
is larger as decreases toward C . g(r) at C is computed by 
summing over all trajectories for all > C . 



We do not expect any singularity of f(u) in u = 16 , and therefore: 



6.4 Isostaticity, g(r) and thermodynamics 

In our argument above the divergence in g(r) does not appear through a 
minimization of a free energy or volume, but rather as a consequence of 
the specific dynamics that took place to reach the jamming threshold. Em- 
pirically, it turns out that the divergence in g(r) was observed in systems 
obtained with different dynamics. It was first noticed in molecular dynamics 
simulations of elastic spheres with friction, where the jammed states were 

16 We do not exclude that there is a finite probability po for particles that lost their 
contact at a given (f>o to recover it at a <f>\ such that <p c < <f>\ < <pQ. In our formalism such 
eventuality will not create a divergence in f{u). The incoming flux of contact would be 
compensated by a larger rate of contact opening, and the expression (|fi . 1|) will be corrected 
by a factor -r-^— . 



g{r) 



1 



(6.6) 



(r-l)2 
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obtained using inelastic collisions [53] . It was also recently observed - 
with possibly a slight difference in the exponent — in random close pack- 
ing of hard spheres obtained by increasing continuously the packing fraction 
from the liquid state [JJJ]. In this last example thermodynamic equilibrium 
is maintained until fairly high packing fractions, corresponding to typical 
inter-particle spacing much smaller than the particles radii 17 . On the other 
hand in this hard sphere system the range of validity of Eq. ()6.6|) extends to 
rather large distances, say r — 1 ~ 0.3. Thus, in this system the divergence 
of g(r) at C cannot be the vestige of an interaction that would have taken 
place between pairs of neighbor particles at lower 0, as the system has a fi- 
nite relaxation time and neighbor particles can move arbitrary far from each 
other. Instead, this suggests that the divergence of g(r) obtained in [73] is 
related to thermodynamic properties of a hard sphere liquid near random 
close packing. 

In what follows we discuss why the properties of isostaticity and diver- 
gence of g(r), which were introduced in relation with the mechanical stability 
of the solid phase, may also be connected to the thermodynamic of the liquid 
phase of hard spheres (see also Chapter |HJ). Our argument is in two steps. 
We first argue that the isostaticity and the divergence in g(r) have an inter- 
pretation in terms of density phase space of hard sphere packing. We shall 
show that if an assembly of hard spheres is (i) sub-isostatic, that is with a co- 
ordination smaller than 2d, or (ii) isostatic with a singularity of g(r) in r = 1 
weaker than Eq. ()6.6j) . then it is not a "good" local maximum of density: one 
can build close configurations which are denser. In a second step, we simply 
argue that in the liquid phase, the system prefers to lie in configurations 
close to the denser packing, which is favorable for entropic reasons. Finally 
this suggests that when <p increases toward C , the system becomes isostatic 
and display at least a square root singularity in g(r) at r = 1. 

Consider a sub-isostatic configuration of hard sphere in a box. As we 
discussed in earlier Chapters, this system is not rigid. Therefore, if an in- 
finitesimal pressure is applied at the boundary of the box, the system cannot 
resist to it and yields. Unstable modes collapse until isostaticity is recov- 
ered. Consequently, the volume of the box decreases and one obtains a 
denser configuration of spheres. The same kind of argument can be made if 
one considers an isostatic configuration of hard spheres, again contained in 
a box, that do not display the square root singularity of g(r) in r = 1. Let 
assume temporally that the particles are not hard, but soft, and interact for 
example with harmonic interactions. Such system is again unstable toward 

17 For mono-disperse spheres thermodynamic equilibrium can be obtained below typi- 
cally <p ~ 0.58, which corresponds to an inter-particle spacing h « 0.03 
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P " 




4>1 4>2 * 

Fig. 6.2: Cycle of compression-buckling-decompression starting from a 
packing fraction <pi and ending at a larger packing fraction <p2- 

compression: an affine compression of strain e increases the coordination by 
an amount that scales as 5z ~ J 1 +e g(r)dr. If g(r) does not display the 
square root singularity in r = 1, one finds Sz <C Sz min ~ p? ~ ez. Thus 
unstable modes collapse until the system is enough coordinated, which low- 
ers the pressure. Then if the pressure is brought back to adiabatically, one 
obtains a new isostatic state of higher density, as sketched in Fig. (jd2j) . 

Now we argue that a liquid of hard spheres lies near the dense jammed 
states. More precisely, we define the jammed state corresponding to a liquid 
configuration as the solid obtained by hyper quenching the liquid, that is by 
increasing rapidly the radii of the particles until a solid is formed. We argue 
that the densest jammed state are favored, as they offer a larger free volume 
per particle in the liquid configuration. Consider for example two jammed 
states of packing fraction <p a < <p b , and a liquid at packing fraction < a . 
If the liquid lies near the jammed state a the free volume per particles varies 
as (j> a — <p. Thus if p a is the probability for the liquid to have a as the jammed 
state, one can estimate p a /pb ~ (^E^) N ■ Thus the densest jammed states 
are favored as the packing fraction increases. For a quantitative treatment of 
the liquid equilibrium, one should enumerate all the possible jammed states 
that are accessible to the dynamics — as there are configurations such as 
the crystal which are not visited in a reasonable time — and compute their 
structure. A similar approach was followed recently in [80J. Our qualita- 
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tive argument suggests that the isostatic states with the divergence in g(r) 
evoked above are good candidates to dominate the dynamics near the glass 
transition, since they are (i) good (in the sense stated above) local maxima of 
density and (ii) plausibly numerous, as we exhibited a simple way to create 
them. If so, g(r) in a random close packing obtained from the liquid phase 
must display a square root divergence. 



7. Elastic response near the jamming 
threshold 



In this Chapter we first derive a few properties of the response to a local 
strain, such as the deformation of one contact, near the jamming threshold. 
We show that at the isostatic limit, such a response spreads out uniformly 
in the entire system. This implies that the soft modes of an isostatic state 
obtained by cutting one contact are extended, as assumed in Chapter El 
When the coordination increases, we show that the energy cost induced by 
such a local strain grow as 5z. 

This suggests that the shear modulus also scales as 5z, as it is observed 
numerically. We confirm this behavior by computing the response of the 
system to a global compression and a global shear. An assembly of repulsive 
short-range particles near jamming behaves as a gel in the following sense: 
the shear modulus is much smaller than the bulk modulus, as it was observed 
numerically (HHl EJ] • In a gel, the response to compression is the one of liquid 
of monomers, whereas the response to a shear is related to the stretching of 
the polymer chains. The cause for the discrepancy between these two moduli 
is different near jamming. If the responses to compression and shear were 
purely affine, then the two moduli would scale in the same way. As we 
shall show the bulk modulus is simply what one would expect from an affine 
response of the system. In the case of a shear, the non-affine displacements 
lower the energy dramatically, which even vanishes at the transition. This is 
possible because the low-coordination allows particles to rearrange without 
much energy cost. Thus this system is a useful tool to study non-affine 
displacements, that might play an important role to describe material failure 
|82j, and that were observed in various systems such as granular matter or 
biological materials, and in simulations of glasses and foams, see e.g. [S3] 
for references. We discuss very tentatively the possible length scales that 
characterize such displacements. 
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7.1 Formalism 

In this section we study the relation between forces and soft modes. We show 
that the equation of force balance is the dual of the equation that defines the 
soft modes. This will be used in the next section to compute the response to 
a strain. 

7.1.1 Force propagation 

Many of the properties discussed below concern the response of the system to 
external forces. It proves convenient to consider our system under the influ- 
ence of an arbitrary set of forces Fi acting on all particles i. At equilibrium 
the sum of the forces on each particle % is null: 



where is the compression in the contact < ij >, the sum is on the particles 
in contact with particle i, and Fi the external force applied oni. In term of 
sign convention we recall that is the unit vector going from i to j. This 
linear equation can be written: 



|f) is the vector of contact tensions and has N c components. F is the vector 
of external forces. Its dimension is (Nd — d(d + l)/2). Indeed there are 
d degrees of freedom for the external force on each particle, which brings 
the term Nd. The term d(d + l)/2 corresponds to the constraints on the 
total torques and forces that must be zero at equilibrium (in what follows 
our notations shall be lower-case for the contact space of dimension N c , and 
upper-case for the particles positions space of dimension Nd — d(d + l)/2). 
Therefore T is an Nd — d(d +l)/2xjV c matrix. 

In the case where N c > Nd — d(d + l)/2 equation ()7.2|) is not sufficient 
to fix the contact tensions. One must apply elasticity theory to compute 
the contact forces for a set of external forces. This is done by writing the 
relation between forces and distances between particles in contact (which 
depends on the interaction potential), and by relating the distance among 
particles to their position. This brings the necessary N c — Nd — d(d + l)/2 
extra-constraints. When N c < Nd — d(d+l)/2, there is in general no solution 
for a given force field, and the system yields under a generic external force 
field. For such a system no contact forces f can satisfy ()7.2)1 unless F is 
restricted. In the last case, when N c = Nd — d ( d + l > ; the system is called 




(7.1) 



T|f) = |F) 
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"isostatic", and there are just enough contacts to equilibrate any external 
force field. Note that in this case Eq. (j7.2j) . which does not depend on the 
interaction potential, is sufficient to determine the contact forces. 

7.1.2 Duality between force propagation and soft 
modes 

As we discussed in earlier Chapters, the static equilibrium of the system can 
be expressed in geometric terms, rather than in terms of the forces as in 
Eq. (j7.2)l . The change of distance to the first order between particles for a 
given set of displacements 5Ri is: 

(5Ri — 5Rj).fiij = 5rij (7.3) 

It can be written as: 

S\SR) = \5r) (7.4) 

where 5r = {Sr^} is the set of distance changes for all contacts. If one 
removes the global translations or rotations from the displacement fields, 
which obviously do not change any inter-particle distance, S becomes an 
(Nd — d(d + l)/2) by N c matrix. Its kernel is the space of soft modes we 
introduced in Chapter |21 When iV c < Nd — d(d + l)/2, the displacements 
are undetermined for a given change of inter-particles distances. 

Now we establish the connection among soft modes and forces, which is 
also present in a similar form in j^UEE]. At equilibrium for any displacement, 
the net work done by external forces and the contact forces is zero: this 
assures the stability of the system. Therefore for any acceptable equilibrium 
force field: 

£ 8R, ■ Ft - Sryfij = (5R\F) - (Sr ■ f ) = (7.5) 

i ij 

For a soft mode, br^ = and we are left with: 

(m\F)=Y^8R i -P i = (7.6) 

i 

The soft modes are equivalent to the constraints on the force field one would 
have obtained by solving Eq. (j7.1|) . Each soft mode represents a direction of 
fragility of the system, and the external forces must be orthogonal to them 
to avoid yielding. Furthermore, applying the definitions of S and T in (|7.5|) 
we have: 

(f \S\SR) = (SR\T\f) (7.7) 
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therefore, introducing the transpose notation, we obtain: 

S = T* (7.8) 

7.1.3 Relation with the Dynamical matrix 

In this Chapter we shall neglect the initial stress term in the energy expan- 
sion, as we did in Chapters El and EJ As we discussed in Chapter we expect 
that this term affects only weakly the plane waves. Thus it should lead to 
negligible corrections in the computations of the responses at zero wave vec- 
tor such as a shear or a compression. In the harmonic case we discuss, we 
have then: 

SE= £ \5rl = hs5R\S5R) (7.9) 

<ij> 

Therefore one finds for the dynamical matrix defined in Eq. (j2.2j) : 

M=S f S (7.10) 

7.2 Relation between the response to a 
strain and forces 

The response to an imposed strain is given by the minimization of the energy 
with respect to the displacement fields. In this section, using simple linear 
algebra and the duality between force and soft modes studied above, we show 
how this minimization can be written as a sum over the vectors in the contact 
force space that satisfy force balance. This will allow in the next sections to 
(i) derive exact results of the response to a local strain (ii) derive the bulk and 
shear moduli assuming well-known empirical fact of force properties, that we 
recall at the end of this section. 

As we discussed in earlier Chapters, in the approximation where the initial 
stress is neglected, the system is equivalent to a set of point particles inter- 
acting with springs. In order to study the elastic behavior of such system, 
it turns out to be convenient to consider the responses that follow arbitrary 
changes of rest length of these springs. This is in fact equivalent to imposing 
dipoles of force. Consider for example two particles i and j in contact, and 
increase the rest length of their spring by an infinitesimal amount j/y. It 
is equivalent to impose a dipole of force where opposite external forces are 
imposed on i and j with Fj = — Fi = UijUij. As we shall see, the response to 
shear or compression can also be easily expressed in terms of changes of rest 
length. 
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We impose an infinitesimal change of rest length on every contact y = 
{i/ij}. Following Eq. (j7.9j) the energy and the displacement field are given by 
the minimization of: 

5E = -mm(S5R-y\S5R-y) (7.11) 

Obviously if S was spanning its image space, we would have 5E = : one 
could always find a displacement |<5R) that leads to a change of distances 
between particles in contact exactly equal to |y). As we said, S is a iV c x 
Nd - d(d + l)/2 matrix. If N c < Nd - d(d + l)/2, S indeed spans its 
image space, and the energy associated with any strain |y) is zero: the 
system is floppy. In the other case, if N c > Nd — d(d + l)/2, there are 
N c — Nd — d(d + l)/2 = Nbzjl relations of dependency among the columns 
of S. One can choose a basis of NSz/2 vectors |a p ), with 1 < p < N5z/2, in 
the space of \8r) such that: 

(a p \S = (7.12) 

The |a p ) are orthogonal to all the vectors <S|5R), for any displacement field 
|<5R). Transposing this relation we have: 

T|a p ) = (7.13) 

which indicates that all the vectors in the space of the |a p ) satisfy force 
balance without external force (|7.1|) . but no others. The |a p ) live in the 
contact-force space, and henceforth we shall denote them |a p ) = |P) = {ffj}- 
In the following we consider an orthogonal unit basis: 

(pip'> = x;^ = v ( 7 - 14 ) 



We can decompose any |y) as: 



y ± + J2 ( fP ly) fP ( 7 - 15 ) 



p=l...N% 



y x , the part of |y) orthogonal to the |f p ), is spanned by the matrix S, 
and therefore does not contribute to the energy when the minimization of 
Eq. ([7.1ip is performed. In other words, there is a displacement field which 
leads to a strain y- 1 . The energy that results of the minimization of Eq. (j7.11j) 
is then simply given by the distance square between |y) and the image vector 
space of S, that is ||y — y^W 2 or: 

&E = \ £ (P|y) 2 (7.16) 
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Furthermore the response to such strain condensates an energy SEij in the 
contact < ij > that follows: 

SE ij = km - yjj? = [ E (f p ly)/^] 2 (7-17) 

Expressions (17.16)1 and ()7.17j) furnish explicit solutions for the minimization 
of Eq. (|7.11jl . They do not give access to the displacement field 5R that follows 
a given strain (that can be computed in principle by inverting S), but to the 
energy condensed globally or in any contact. From these expressions we shall 
derive a few exact properties of the response to a localized strain or force 
dipole. Furthermore, they relate the response to a strain to properties of 
forces, which are well studied empirically. Making simple assumptions on 
these contact force fields, we shall derive the scaling of the elastic moduli. 



7.3 Response to a local perturbation 

We impose a local strain corresponding to the compression or the stretching 
of one contact. We implement it by increasing the separation r 12 of a single 
contact 12 chosen to be at the origin by an amount e. This is equivalent to 
impose a force dipole Fi = — F 2 = eii\i. We have: 

yij = if ij ^ 12 (7.18) 
2/i2 = e (7.19) 

It is now straightforward to compute the energy cost of local strain using 
Eq.fEEJ. One finds: 

5E = \ £ {f[ 2 ef (7.20) 



1=1,.. N^f 



If one sums this expression for every possible contacts ij and use the nor- 
malization of the vectors of force contact, one finds exactly N^fe 2 . Therefore 
we have the following exact result for the average energy induced by such 
deformation: 

(SE) = ge 2 (7.21) 

This results can be extended to the case where several contacts are deformed. 
For example, if a particle swells, all its contacts are compressed by the same 
factor e, and one finds that the energy still goes as 5z. In a isotropic continu- 
ous elastic medium with a weak shear modulus G, the energy resulting from 
the swelling of a small sphere goes as G. Thus we may infer from Eq. ()7.21j) 
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that G ~ 5z. This scaling is observed numerically [5U]. In the next section 
we shall justify further this scaling law by computing the response to a global 
shear. 

Eq. (j7.21j) shows that when the jamming threshold is approached, deform- 
ing a contact (or imposing a force dipole) becomes softer and softer. Even- 
tually, when 5z — > 0, the deformation of the contact becomes totally soft. 
The response corresponds then exactly to the soft modes that we described 
in Chapter El that appear when a contact (here the contact 12) is removed 
in an isostatic system. A fundamental question was the extension of such 
modes, that was assumed to be over the whole system at that point. In what 
follows we show that this is indeed the case. 

When the jamming threshold is approached, the coordination diminishes 
until there is eventually only one term left in the sum of Eq. (j7.2()j) . At that 
point there is one more contact than in an isostatic state 18 and Sz — jr. 
Eq. (j7.21|) indicates that the energy resulting from the deformation of one 
contact is of order 1/N, which vanishes as the system size increases. The 
deformation becomes soft. We can use Eq. (j7.17|) . whose sum contains only 
one force field, to obtain the spatial repartition of the energy: 



jf 1 } is the physical set of contact forces that support the system. These 
contact forces are well studied numerically and experimentally, and we shall 
discuss more their properties in the next paragraph. At large distance where 
the contact force are de-correlated ((f^flj) 2 ) = ((fij) 2 ) 2 > 0: the energy 
condensed in the contact ij does not vanish when the distance between ij 
and 12 diverges. If the displacements that follow such a perturbation were 
to vanish with the distance from the source, it would be so for the energies 
condensed in the contacts. As it is not the case, these displacements spread 
out in the entire system. 

18 At the jamming threshold, the system is in fact not exactly isostatic: it has one extra 
contact. It must be the case as there must be a non-zero contact force field to pin the 
particles. The counting of degrees of freedom of Chapter N c — Nd — d(d + l)/2, is in 
fact incomplete: it must take into account the degree of freedom corresponding to the box 
size. This subtle point, which does not affect the result of the previous Chapters, is trivial 
in one dimension, as N particles pinned on a ring have N contacts, and not N — 1. This 
was verified numerically in the simulations of |50j in 2 and 3 dimensions as removing m 
contacts at the threshold creates only m — 1 soft modes instead of m. 
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7.4 Elastic moduli 

7.4.1 Spatial properties of the force fields \f p ) 

In this paragraph we discuss the properties of the contact force fields that 
we shall use to derive the scaling of the bulk and the shear moduli. Only one 
vector of the vector space of the force fields |f p ) solutions of Eq. (j7.2j) without 
external force is the real set of contact forces that supports the system. As we 
discussed at the beginning of this Chapter, it can be computed by solving the 
whole elastic problem. The solution of such problem is unique. This vector 
is denoted jf 1 ). The rest of the basis |f p ) with p ^ 1 are also solutions of Eq. 
()7.2|) without external force. Nevertheless, there are not "physical" solutions 
for the interaction potential chosen. Thus we shall call them virtual. 

|f x ) verifies the following properties: (i) In a system with repulsive in- 
teraction, as we consider here, all the contact forces are compressive and 
therefore > for all contacts, (ii) It is well known from simulations and 
experiments that the distribution of contact forces is roughly exponential, or 
compressed exponential (see for example [SUj for simulations in the friction- 
less case). This implies that the fluctuations of the contact forces are of order 
of the average value, leading to (flj) 2 ~ ((flj) 2 ) — 1/N C for a normalized force 
field. Thus we may introduce a constant Co such that: 



Now we turn to the properties of the virtual forces |f p ): (i) There are no 
physical constraints on the sign of the contacts forces for the virtual vectors. 
Furthermore, the |f p ) must be orthogonal to |f 1 ), whose signs of contact forces 
are strictly positive, and where the fluctuations in the contact compression 
is small. Therefore the virtual force fields have roughly as many compressive 
as tensile contacts. 

7.4.2 Implementation of global strain 

In our framework it turns out to be convenient to study the response to shear 
or compression as there are generally implemented in simulations. When 
periodic boundary conditions are used, an affine strain is first imposed on 
the system. Then the particles are let to relax. In general the affine strain 
is obtained by changing the boundary condition. Consider a 2-dimensional 
system with periodic boundary: it is a torus. For example a shear strain 
can be implemented by increasing one of the principal radii of the torus 
and decreasing the other. Then the distance between particles in contact 




1 



(7.23) 
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increases or decreases depending on the direction of the contact. In fact, 
this procedure of change of boundary conditions is formally equivalent to a 
change of the metric of the system. If the metric is changed from identity 
I to the constant metric G = I + U, the length of a vector 5lo becomes SI, 
such that SI 2 = SIq ■ G ■ 5l . From this expression one deduces the change of 
distance between two particles is given by the formula: 

5 rij = Hij ■ U ■ Rij (7.24) 

Near jamming R^j w and therefore Sr^ w fty ■ U ■ n^. Formally, such a 
change of metric is strictly equivalent to a change of the rest length of the 
interactions with = fiy • U ■ n^. Incidentally Eq. (|7.1fij) can be used to 
compute the energy of such strain. 



7.4.3 Compression 

For a compression U = —el where I is the identity matrix and e is the 
magnitude of the strain. Eq. (|7.1fij ) becomes: 

8E = lCE-e%)* + l £ (-eT,& (7-25) 

In the first sum all the terms have the same signs for a purely repulsive 
system, and this term leads to the strongest contribution. We have: 

SE > (E "^) 2 = t\N c {k)f = e 2 c 2 N c (7.26) 
u 

On the other hand, 5E is certainly smaller than an affine compression 
whose energy also goes as e 2 N. Therefore we find that: 

5E ~ A^e 2 (7.27) 
B = ^~6* (7.28) 

The bulk modulus of an harmonic system jumps from in the "gas" phase 
toward a constant when the system becomes jammed, as verified in the sim- 
ulations. From this result follows that p ~ (<f) — C ). Note that this result 
holds only for purely repulsive systems. If there are as many tensile and 
compressive contacts, one recovers for the bulk modulus the result valid for 
the shear modulus, that we derive in the next section. 
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7.4.4 Shear 

If a pure a shear strain is imposed, the tensor U is traceless. Let be e the 
largest eigenvalue (in absolute value). The change of distance of two particles 
in contact due to shear br^ = ■ U ■ n^. It is a number of zero average if 
the system is isotropic, and fluctuates between +e and — e depending on the 
orientation of Rij. Eq. (j7.16|) becomes: 

*E = \ £ (7-29) 



2 



11 

2 J 



Each term in the summation gives on average: 

((£ Jtfir*?) = E((& S 4) + £ UlfL»r ; ,6r ) (7.30) 

ij ij mnj^ij 

= £<(/£) 2 ><<H>+ E m P mJra5r mn ) (7.31) 

In principle they can be spatial correlations in the contact force ampli- 
tudes that leads to (fijfmn^ r ij^ r mn) 7^ even if mn ^ ij. Nevertheless we 
expect these terms to be negligible, as we argue at the end of the paragraph. 
Concerning the diagonal terms, one has Srfj fa e 2 while J2{fij) 2 = 1 by 
construction. Thus each term in the p summation is of order e 2 • 1, and: 

5E ~ N5ze 2 (7.32) 

which implies: 

G=^~6z (7.33) 

This is in agreement with the observation of [50] . 

We come back to the subtle question of the possible presence of spatial 
force correlations. In the simulation of soft spheres of [HO], the spatial corre- 
lations of the real force field |f x ) vanish for distance larger than roughly two 
particle diameters at the jamming transition. This property is restituted by 
simple models of force propagation such as the scalar g-model [HI]. In this 
model, the system is decomposed in layers, and forces propagate downward 
from the top of the system to the bottom. Each particle as 2d contacts, 
with d particles of the upper layer and d particles of the layer below. The 
force field is builded recursively, following a local rule: when a particle re- 
ceives from an upper layer a total force amplitude / = fi + ■■■ + fd, this 
force is distributed randomly to the d contacts below. This model mimics 
the fact that in an isostatic system external forces imposed at the top of a 
system propagate downward, as we discussed in Chapter El The randomness 
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of the force splitting sketches the randomness of the local configuration of 
contacts. This model presumably captures some of the relevant physics, as 
it leads both to a rather realistic force distribution, and does not display any 
spatial correlation in the force amplitude. Extending this model, we argue 
that when 5z > the typical virtual force fields do not display long-range 
correlations. Indeed following the same line of thought a typical virtual force 
fields can also be builded recursively with a similar local rule, the only differ- 
ence being that some particles have more than 2d contacts. Thus the level 
of randomness increases as there are sometimes more ways to split the force 
between the different contacts below, accounting for the fact that there are 
many possible virtual force fields. Such an increase of randomness will cer- 
tainly not create long-range correlation, and this simple model justifies our 
assumption that such correlations are negligible. Note that this model does 
not preclude that a few force fields can display long-range correlations, but 
simply supports that such correlations do not occur for the bulk part of the 
set of force fields. In particular, we expect that when 5z > 0, the real force 
field If 1 ) displays long range correlations. Indeed this force field is solution 
of the whole elastic problem, and cannot be builded from any local rule. We 
expect that at large distances r, the correlation in force amplitude would 
follow (/ 1 (0)/ 2 (r)) ~ r~ d as in a continuous elastic medium with random 
disorder jHH|. Accounting for this particular force field in Eq. (J7.30J) leads to 
a relative correction of G that goes as ln ffl which vanishes at the thermody- 
namic limit. Finally, note that some subset of these virtual force fields were 
studied in simulations [S3] , and that no long range correlations were noticed 

159 - 

7.5 Discussion: non-affine displacements 
and length scales. 

An affine shear costs an energy comparable to a compression. Thus the non- 
affine displacements that follow a pure shear diminish the cost in energy by a 
ratio that diverges at the jamming transition. Such non-affine displacements 
I^R-n.a.) are simply the displacements corresponding to the minimization of 
Eq.fJTTTJ). It follows that ly -1 ) = <S|<5R n . a .). S can be made invertible if it is 
restricted to its image space, and we may write: 

\6R n . a .) = |5-V> (7-34) 

A surprising recent observation was made in the Lennard- Jones simulations 
of jHllZ|: a rather large length scale (~ 30 particle sizes) appears in the cor- 
relations of the non-affine displacements. Eq. ()7.34j) is a linear equation. Any 
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strain |y) can be decomposed as the sum of individual contact deformations. 
Therefore the non-linear displacements that follow a shear or a compression 
can be written as the sum of the responses to a contact deformation. The 
formalism of the present Chapter does not give any access to the spatial be- 
havior of the response to a local deformation when 5z > . On the other 
hand the results on the vibrational modes of Chapter 0] might bring interest- 
ing insights. Since any deformation can be decomposed on the vibrational 
modes, we expect the characteristic lengths I* and k to characterize the non- 
affine displacements that follow a global strain. It would be useful to study 
these questions numerically near the jamming threshold. 
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In the previous Chapters we discussed some properties of weakly-compressed 
harmonic soft spheres. Our main results are that (i) anomalous modes appear 
at low-frequency. They are related to the soft modes of floppy subsystems 
and are characterized by a length scale I* (ii) there is a frequency scale u* 
below which the system does not behave as a continuous medium, but as 
in an isostatic state. In this Chapter we study the applicability of these 
results to more realistic systems. We first discuss granular matter, then 
glasses. The sections contained in this Chapter, apart the last one, can be 
read independently since they deal with distinct matters. 

At least two modifications are necessary to apply our results to granular 
matter: (i) the presence of friction. Up to now we studied mostly radial 
interactions, and covalent bonds in Chapters El and El In section 18.11 we 
shall see that there is no conceptual difference when friction is present, apart 
from a change in the equations that define the soft modes, (ii) the potential 
used: In three dimensions grains do not interact harmonically. Modeling the 
contact with a Hertzian potential is more realistic, even thought it might 
not be perfect [S*7| . It corresponds to a = 5/2 in Eq. (|l.lj) . In section EZ1 
we compute the density of states for such non-harmonic interactions. Our 
results agrees with the numerical findings of |5D] . 

Then we discuss the vibrations of glasses. All glasses present attractive 
forces, such as Van der Waals interactions. Therefore each particle a pri- 
ori interacts with all other ones, and the coordination number is not-well 
defined. In section 18.31 we show how one can deal with this problem and 
compute the density of states in simple cases where there is a strong hierar- 
chy in the contact stiffness. Our results apply for systems as silica, where the 
covalent interactions inside the tetrahedron are the strongest as discussed in 
Chapter 01 In section 18.41 we consider systems without any clear-cut hierar- 
chy in the distribution of contact stiffness, such as Lennard- Jones systems. 
The question is to know under which conditions our results on the vibra- 
tions of weakly-connected system can apply in such situations. We propose 
an improved variational argument which uses the distribution of the contact 
stiffnesses to evaluate the density of states. This leads to testable predic- 
tions on the nature of the excess-modes in such systems. In particular, this 
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should enable to decide whether or not the length scales that appears in the 
vibrations, the responses to a point force or the non-affine displacements in 
some Lennard- Jones simulations [HI IE1 EI] corresponds to the length scales 
we introduced in this thesis. 

Let us specify that in this Chapter we shall not consider the effect of 
initial stress: this is a separated issue, which can be treated independently, 
as shown in Chapter 

8.1 Particles with friction 

In this section we show how our results can be extended to particles with 
friction. We shall assume that all the contacts lie inside the Coulomb cone: 
particles in contact do not slide irreversibly for an infinitesimal perturbation. 
This means that there is no plastic deformation at the surface of contact: the 
contacts act as if they were "welded". Simulations showed that the validity 
of this assumption depends on the system preparation. The distribution 
of the contact tangential force was observed to vanish on the edges of the 
Coulomb cone in [HHI- In [88J, this was also observed if the dynamics that 
leads to jamming is fast. Otherwise, a finite fraction of the contacts was 
observed to lie on the Coulomb cone. For concreteness we consider the two 
dimensional case where the notations are simpler. The same ideas apply in 
higher dimensions. We consider elastic discs, whose repulsive interaction is 
harmonic: a = 2. We add a term of friction, function of the shearing of 
the contacts. If two particles are adjacent, and one of them rotates an angle 
SO while the other is rigidly pinned, the energy stored can be written in the 
form: 

5E = 1 56 2 (8.1) 

7 describe the energy associated with the shearing of a contact. This can be 
easily generalized in 3-dimensions [89J. 

If a small displacement |5R) and a small rotation field 56 = {56i} are 
imposed, the shear between two particles in contact is characterized by the 
following angle: 

56^ = 6 { + Q 3 + (6Rj - SRj.nt- (8.2) 

where njj is obtained by rotating by +|. Therefore, using (j8.1j) . we can 
write for the expansion of the energy: 

5E = E \[( 5 Rj ~ 5 Ri)-nij} 2 + 7ft + Oj + (SRj - SRi).^} 2 + 0(5R 3 ) (8.3) 

ij 
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For each contact ij there is a compression force and a tangential force 
fjj. The system is at equilibrium when the force and the momentum are 
balanced on every particle: 

E^i + /iH = ^ foralli ( 8 - 4 ) 

J2fi = Mi foralli (8.5) 

j 

where Aii is the external momentum applied on particle i. This set of linear 
equation has 2N C degrees of freedom and 3N — 3 constraints in two dimen- 
sions. Therefore a jammed system with friction that can sustain any generic 
external force field must be such that 2N C > 3N — 3. At the isostatic point 
the coordination number is given by: 

z c = 2NJN -> 3 (8.6) 

In three dimensions one finds similarly z c = 4. The soft modes are modes of 
null energy. In two dimensions Eq. (j8.3)l gives: 

{5Rj-5Ri).n i3 = U (8.7) 
Oi + 9 j + (SRj - 6Ri).n± = (8.8) 

Again this is a linear system. Now there are 3iV — 3 degrees of freedom, and 
2N C constraints. Following the procedure of Chapter El and |U one can build 
anomalous modes that appear at a frequency uj* ~ 5z, where 5z = z — z c , 
z c being the isostatic limit with friction. This result was recently observed 
numerically in [HOI- Furthermore, it is easy to show that this system (J8.7)) is, 
as in the frictionless case, the dual of ([8.4)1 . The results on the dipole of force 
and the scaling of the elastic moduli of Chapter [7| can then be recovered. 

Note that contrarily to the frictionless case, this system does not need 
to be isostatic when the pressure vanishes. As we discussed in Chapter El 
and El in frictionless systems the two conditions that (i) the system must 
be rigid and (ii) the particles cannot interpenetrate lead to the unique solu- 
tion for the coordination number: z = z c . In the frictional case, these two 
antagonist requirements do not lead to a unique solution. The geometrical 
requirement that forbids the interpenetration of particles does not change, 
whereas the condition of rigidity becomes less demanding, as we just showed. 
In two dimensions one finds 3 < z < 4, and in three dimensions 4 < z < 6. 
Consequently the coordination of stiff frictional particles depends on the 
preparation of the system [SHI ESI- I n three dimensions when the friction co- 
efficient [L m 1 the simulations of [HBJ leads to z ~ 4.5 that implies Sz ~ 0.5. 
We guess that this corresponds to a /* of few tens of particle sizes, which 
could be probed by computing the vibrational modes of this system. 
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8.2 Extension to non-harmonic contacts 

Up to now we considered harmonic interactions. Here we discuss the general- 
ization of our argument to other contact potentials. Ref. [SU] explored several 
other potentials, in particular the Hertzian contact potential describing the 
compressive energy of two elastic bodies. It corresponds to a = 5/2 in Eq. 
(jl.lj) . Ref.jnm observed a plateau in the density of states whose height D 
scales as p~ 1//6 . They also observed a cutoff frequency uj* varying as pz. In 
the Hertzian case the quadratic energy of Eq. (|2.1|) becomes: 

*E=\ E(l - ^[(5% - *&).fjy] a (8.9) 

The new factor (1 — r^)^ amounts to a spring constant fey that depends on 
compression. The contact force = dbEjdr^ evidently varies as (1— r^) 3 ' 2 . 
In what follows we start by neglecting the fluctuations that exist between the 
stiffnesses of the contacts. This treatment is sufficient to recover the scaling 
results of [201 ■ Then we discuss how such fluctuations can be taken into 
account to improve the bound on the density of states derived in Chapter El 
The new factor V"(ry) = (1 — r^)^ rescales the energy. To account for 
this overall effect, we replace (1— by its average ((1— r^)^}. Expressed in 
terms of contact forces, this factor is proportional to (fij 3 )- Again replacing 
fij by its average, the factor becomes (fij) 1 / 3 - This average is related to the 
pressure p, via p ~ (fij)- Thus in this approximation the overall effect is to 
rescale the energy by a factor k(p) ~ p 1 ^ 3 . 

Apart from this prefactor, the energy and the dynamical matrix are iden- 
tical to the harmonic case treated above. Thus each normal mode frequency 
gains a factor k^ ~ p 1 ^ 6 . In the harmonic case the crossover frequency 
follows uj* ~ 5z. In the Hertzian case, it too gains a factor fea ; so that 
uo* ~ k?5z. When the initial stress is taken into account, the bound on the 
lowest-frequency anomalous modes u>am still has the form 

ujam 2 < oo* 2 - A 2 p (8.11) 

For a marginally-stable system we still have u>am ~ which leads to an 
unaltered relationship between uj* and p, namely uo* ~ p%. Comparing with 
our previous estimate of uj* yields 5z ~ p 1 ^ 3 . Furthermore, the plateau 
density of states Dq has the dimensions an inverse frequency, and thus gains 
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a factor p -1 / 6 . Since the harmonic D had no dependence on p, the Hertzian 
Do(p) also should vary as p -1 / 6 . The scaling behaviors seen in jHOj agree with 
these expectations. These arguments may be applied to general values of the 
interaction exponent a. 

Additional effects could in principle alter the low-frequency modes in the 
Hertzian case. When harmonic springs are replaced by Hertzian springs, the 
contacts supporting different forces have different stiffnesses. The variational 
argument derived in Chapters El and |U does not consider such fluctuations. 
If it is applied as is to the more general case with fluctuations, the energy of 
the corresponding anomalous modes defined in Eq. (|3.1|) simply gains a factor 
(k), the average stiffness. This can be deduced from Eq. ()3.3|) by neglecting 
the correlations between the soft modes displacements and the stiffnesses 
amplitudes 19 . Thus, this variational argument corresponds to the simple 
derivation of the previous paragraph where quantities are replaced by their 
average. Here we propose to use the stiffness fluctuations to improve this 
variational argument by a numerical factor. Instead of modulating the soft 
modes by sin(^) to obtain the anomalous modes of Eq. (j3.1|) . we introduce 
a more general phase and modulate the soft modes with sin(^(z)). Then 
we minimize the average energy of these anomalous modes with respect to 
ip, imposing ip = and ip = ir on the boundaries x = and x = L. The 
expression of the energy corresponds to Eq. ()3.3|) with sin(^>(z)) replacing 
sin (^jr-) and each term of the sum multiplied by a stiffness fey. When averaged 
on the soft modes amplitude, one obtains the average energy of the anomalous 
modes (5E) ~ YUj hj cos 2 (ip(i))[ip(i) — ip(j)] 2 . This has to be minimized with 
respect to the variables We expect the phase tp to vary slowly in space, 
and not to differ dramatically from our previous solution <fi = ttx/ L. Hence 
for a large system we may consider that the cosines terms that appear in 
the energy sum are constant locally. Then the local minimization of the 
energy corresponds to the problem of conductivity in a random network 
of resistors where SE = Y.ij r ij~ 1 [U{i) — U(j)} 2 , where is the resistor 
between the vertex i and j, and U(i) is the Coulomb potential in i. Therefore 
kij corresponds to the inverse of a resistance r^, and ip(i) to the potential 
U(i). Effective medium theory j^I] furnishes a good approximation of the 
conductivity, and therefore of the energy 5E, of such a random system. For 
example if the network of contacts is sketched by a square or cubic lattice, 
one obtains for the effective stiffness the following equation: 

< (d -XW =o (8 - 12) 



19 We expect such correlations to be very small, as the amplitudes of the stiffnesses do 
not enter in the soft modes equation. 
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where the average is taken on the distribution of stiffness. Because of the 
convexity of the inverse function, one obtains that fc e // < (A;). Expectedly, 
the energy in this improved variational argument is smaller than in our previ- 
ous result, since it is proportional to fc e // instead of (k). For the distributions 
of stiffnesses we are considering here 20 where we expect no delta function in 
k = 0, nor any fat tail at large k, k e ff as computed in Eq. (j8.12|) and (k) are 
of the same order of magnitude. Thus the bound on the density of states is 
improved by a numerical factor when this improved argument is used. 

8.3 D(uj) in systems with various interaction 
types 

In this section we study the density of states of systems where several types of 
interactions are at play, and where the amplitude of the strongest interaction 
is much larger than the others. In particular we think about tetrahedral 
covalent networks such as silica. The strong interaction corresponds to the 
deformation of the tetrahedra, which is much larger than the bending of the 
Si-o-Si bonds or the Van der Waals attractions, as we discussed in Chapter El 
In Chapter|3]and|3]we showed how to compute the density of states of weakly- 
connected networks with only one type of interaction: no weak interactions 
were present. The density of states is then described by a plateau that 
appears above a frequency uj* ~ 8z. In this section we show how to treat the 
weak interactions by simple perturbation. We simply consider the anomalous 
modes that would appear at u* if the weak interactions were not there. Then 
we compute the change in energy induced by the weak interactions on these 
anomalous modes. If these modes are on average shifted by an energy (, the 
plateau in the density of states appears at a frequency yjui* 2 + (. In what 
follows we show how to evaluate ( and discuss some examples. 

8.3.1 Model 

Consider for concreteness a system with radial interactions, where "strong" 
interactions with stiffness close to unity form a rigid network of coordination 
z > z c , and where the weak interactions have stiffnesses of order r] 1. 

20 As we discussed in Chapter the potential does not enter in the computation of the 
force field at the isostatic threshold. Thus if the configurations obtained at threshold are 
similar for the different potentials a, as it seems to be the case, the same force distribution 
is obtained for any a. Such force distribution is known to be roughly exponential. From 

i 

this we can deduce the distribution of stiffnesses P{k) ~ ^r'"" . 
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Consider the lowest-frequency anomalous modes that appear at u* in the 
rigid network of coordination z. The energy correction ( induced by the 
weak interactions on such mode is: 



where the sum is taken on all pairs of particles interacting with the weak 
interaction of stiffness rjij. Here we consider that all pairs of particles < ij > 
that appears in this sum are not in contact through the strong interaction, 
since such terms would simply re-normalize the strong stiffnesses, and can be 
taken into account when uo* is computed in the first place. To estimate (() we 
need to compute ([(5Ri — 5Rj) ■ n^} 2 ). If the displacement of particles i and j 
were uncorrelated, one would get for a normalized mode ([(5Ri— 8Rj)-nij] 2 ) = 
■jjj. We expect this to be the case when the distance between the two 
particles is large. When the distance between the pair ij diminishes, the 
displacement of particles i and j becomes correlated. Thus we may write 
([(8R4 — SRj) ■ riij) 2 } = 2C jx\ where c(ry) < 1 characterizes the correlations 
at distance r^. Since we are only considering pairs of particles which are not 
in contact for the rigid network, we expect that for the anomalous modes 
of such a network considered here, the correlation of the displacement of 
particles i and j is "weak", that is c(rij) is bounded below by a constant 
Co of order one. This is supported by our numerical result of Fig. fj5.lj) . that 
shows that if two displacement degrees of freedom are not fixed by the soft 
mode condition, then the correlation between these displacements is weak, 
even when these two particles are close. At the level of approximation we are 
considering here we shall simply assume that c(ry) = 1. Then we obtain: 



where we introduced the number p{rj) of pairs with stiffness rj per particle 
per unit stiffness. 

8.3.2 Density of states of square and cubic lattices 

In what follows we discuss in particular the situation where the rigid network 
is isostatic z = z c . Our present argument shows that a plateau appears in the 
density of states at u> ~ To illustrate this idea we consider the example 
of the square (or cubic) lattice of point particles where first neighbors interact 
with spring at rest, see Fig. fj8.lj) . Our argument of Chapter which does 
not depend on disorder but only on coordination, applies to this system. 




(8.13) 
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The square and the cubic lattice are isostatic, therefore there density of 
states is constant. In this particular case the soft modes are trivial one- 
dimensional object: they are the independent translations of the columns. 
Thus for these systems a simpler way to prove that the density of states is 
constant, instead of going through the argument of Chapter El is to consider 
such lattices as an assembly of disconnected one dimensional lattices, which 
obviously have a constant density of states. One may ask what happens to the 
density of states of such system if a weak interaction is added to couple these 
independent lines. As an example, we consider that the second neighbors are 
now coupled: we add springs with small stiffness rj as drawn in Fig. (j8.1|) . p is 
then a delta function at stiffness 77, of amplitude 2 if one considers the square 
lattice (there are two weak interactions per particles). Therefore our rough 
evaluation for the appearance of the plateau of Eq. ()8.14|) yields uj = ^Jf\. 
Since this system is crystalline, the normal modes are plane waves and the 
density of states can be computed exactly. The density of states display a 
cross over between a Debye regime and a plateau. A direct estimation of the 
cross-over energy gives 5E « 2rj which leads to u — \fh] 21 . A cross-over 
energy is arbitrarily defined, here we simply note that our two estimates are 
similar. Finally, note that our argument indicates that the cubic lattice are 
very floppy if the second neighbor interactions are small in comparison with 
those of the first neighbors. This happens if the interactions potential decays 
rapidly. In practice, there are no simple elements that crystallize in a simple 
cubic crystal j2Hj because such structures are too floppy. For example it is 
not possible to crystallize a 12-6 Lennard-Jonnes in a square lattice even at 
zero pressure and zero temperature, because such a structure is mechanically 
unstable 22 . In general simple cubic crystals have charged particles, as in 
NaCl, and the second neighbors interactions are not negligible. 

21 The frequency of a plane wave of wave vector k can be written u> ~ cl(k/k) sin(«;), 
where the function a(n/ k) only depends on the direction of the wave vector. If the sine is 
approximated by an affine function on obtains to ~ a(H/ 'k)k. Then it is straightforward to 
show, by summing on all directions k/k, that in this approximation the density of states is 
linear (as expected for a 2-dimcnsional crystal) until a frequency too, defined as the lowest 
frequency at which some plane waves reach the edge of the Brillouin zone. These plane 
waves have the smallest a(n/ n). In this example they correspond to the transverse waves 
whose wave vector is in the direction of the axis x or y as indicated in Fig. (|8. lfl . When 
the wave vector reach the boundary of the Brillouin zone, the energy of such waves can 
be computed. One finds E — 2r/. 

22 Imposing that the pressure is zero brings a condition between the forces carried by 
the first neighbors and by the second neighbors. In a 12-6 Lennard-Jonnes it can be sat- 
isfied only if the second neighbors are located behind the inflection point of the potential. 
Therefore the second neighbors have a negative stiffness, which destabilizes the marginally 
stable network constituted by the first neighbors. 
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Fig. 8.1: Left: Square lattice with interaction between first neighbor of stiff- 
ness 1. Right: the same lattice where interactions of stiffness 
T] ^ 1 have been added between second neighbors. 

8.3.3 The boson peak of silica 

Silica is the most common glass, much studied in experiments and in simu- 
lations. It is also known to have one of the strongest excess of low-frequency 
modes, or boson peak, see jlHj for a review of empirical results and mod- 
els. In this paragraph we propose an explanation for its density of states 
at low-frequency, and for the nature of the excess-modes. As we discussed 
in Chapter the strongest interactions in silica lies inside the tetrahedron 
SiCv If the other interactions are neglected, one can model the glass as an 
assembly of linked tetrahedra of appropriate stiffness loosely connected at 
corners, the "rigid unit modes" model [10] • This description is supported 
by recent Hyper-Raman scattering experimental results that show that the 
low-frequency excess modes correspond to the motion and rotation of stiff 
tetrahedra [53]. As we discussed in Chapter H] such a tetrahedral network has 
a constant density of states at low frequency. Following the results of the last 
section, one may evaluate the effect of the weaker interactions (in particular 
the bending of the Si-o-Si bond and the Van der Waals interaction) on the 
density of states by computing the energy shift ( that they induce on the 
anomalous modes. ( can in principle be estimated from the knowledge of the 
distribution of the stiffnesses of the weak interactions p. Using the stiffness 
of the Si-O-Si bending interaction obtained ab initio [12], and the molecular 
mass, we obtain a frequency 1.4 Thz, which estimates the order of magnitude 
of a/C. According to the previous arguments of the present section, we ex- 
pect thus silica glass to display a plateau in its density of states that should 
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appear at a frequency of the order of 1 Thz. This is indeed what is observed 
in simulations: silica glass present a well-defined plateau in the density of 
states, which appears at a frequency corresponding to the boson peak, see 
e.g. Fig. 15 of [HI] or Fig. (j8.2j) for numerical results. 

Our argument does not involve disorder. Thus it must also apply for 
the crystals of the same composition and similar densities such as the a and 
/3-cristobalite, since these crystalline structures are formed, as silica, by Si0 4 
tetrahedra connected at the corners. /3-cristobalite has the structure of the 
diamond, in which the tetrahedra correspond to the 4 carbons bonded to a 
central carbon, whereas a-cristobalite has a tetragonal structure. Empirically 
a boson peak is observed in all these materials jHSj- Numerically, a plateau 
indeed appears in D(u) at roughly the same frequency in the cristobalite a 
and (3 [HB] and in the glass, as shown in Fig. (|8.2jl . 

More generally, there are crystals showing an excess density of states at 
frequencies that correspond to the typical boson peak frequency in glasses 

Wf\ l4*Hl l4*§] . This implies that the disorder is not a necessary condition 
to obtain excess modes. This is a crucial point as in many theories of the 
boson peak disorder is the only relevant parameter. Such theory certainly 
cannot explain the density of states of silica, very similar to the cristobalites. 
Rather, we argue here that the key feature that determines the density of 
states is weak coordination. 

Note that if disorder is not relevant to compute the density of states, 
it affects the nature of the vibrational modes. For example properties such 
as transport are very different between silica glass, and the cristobalites. 
Thus the peculiarity of the amorphous state lies in the nature of the excess- 
modes, not in the density of states It is useful to note the parallel 
between cristobalite and silica glass on the one hand and cubic lattice and 
the jamming threshold of elastic spheres on the other hand. In both cases 
the amorphous solid and the crystal have a similar density of states, but the 
anomalous modes in the amorphous phase are not plane waves. Disorder 
strongly affects the anomalous modes, and makes them very heterogeneous, 
as it appears in Fie l3.ll 

There is an interesting geometrical parallel between a tetrahedral network 
and an assembly of spherical particles with friction. The equation that defines 
the soft modes of these two systems are the same, and corresponds to the 
extension of Eq. ([8.7|) in three dimensions. This comes from the fact that 
when friction is present, the contacts can be considered as "welded", as 
we discussed in 18.11 Similarly in a tetrahedral network the corners of each 
tetrahedra in contact are attached. In both cases the soft modes corresponds 
to the rigid motion of the particles that keep the contact point welded. 

Finally, we considered systems in which the strong interactions form a 
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Fig. 8.2: Density of states of silica glass (at temperatures of 10 and 300 
Kelvins), a-cristobalite and /3-cristobalite. This figure is taken 
from the simulations of [68] . 
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rigid network, that is 5z > 0. Our finding is that the density of states 
can be understood in terms of the anomalous modes that we introduced in 
Eq. (|H.l|) . Note that if the concept of such anomalous mode is new, soft 
modes have been used in several different fields as we discussed in Chapter 
El When the network of strong interactions is not rigid, 5z < 0, and the stiff 
network has got some soft modes. Consequently there is delta function in 
the density of states at uj = if the weak interactions are neglected. When 
the weak interactions are taken into account, Thorpe and collaborators |HZ| 
argued that these soft modes shift to higher frequency (in our notation by 
an amount \/C)> leading to a peak in the density of states. The existence 
of such peak was observed in neutron scattering data of Ge^Sei-^ glasses. 
When the composition of this glass is changed as x decreases, the covalent 
network becomes less and less connected. At x — 0.24 it is isostatic. The 
peak amplitude was observed to diminish as x decreases toward 0.24 P?Hl97j. 
as predicted. 

8.4 Lennard-Jones systems 

In this section we discuss the vibrations of systems where (i) coordination 
is not well defined a priori and (ii) there is not a clear-cut hierarchy in 
the amplitudes of interactions. We have in mind Lennard-Jones systems, 
or purely repulsive system with a potential of the form V{r) = r _/3 . In 
particular, we discuss under which conditions such systems should display the 
anomalous modes we introduced in Chapters 01 One important motivation is 
the recent results on Lennard-Jones systems 13 El El where excess-modes 
are observed, and where a characteristic length scale appears, in particular 
in the vibrations of the system and in the response to a local force. It is 
tempting to associate the excess-modes of [6 to the anomalous modes, and 
the length scale observed by the authors to the lengths we introduced in 
the previous Chapters. In this section we propose an improved variational 
argument, in the spirit of the last section, to evaluate the density of states 
of such systems. We shall argue that this variational method is efficient 
if the interaction potential decays rapidly. In this case anomalous modes 
appear characterized by a length scale. We propose a numerical test to check 
whether or not the excess-modes found in the Lennard-Jones system of [Oj 
correspond to the anomalous modes we discussed. 

For concreteness we consider a system with radial interactions, and in- 
troduce the density of stiffness P{k), which denotes the number of pairs of 
particles whose interaction has a stiffness k per particle per unit stiffness. In 
the harmonic soft sphere system that we considered in the previous Chap- 
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ters P(k) = ~8(k = 1). In the previous section, P(k) could be decomposed 
in two well-separated distributions, one around k — 1 of amplitude z/2 for 
the strong interactions, and one at much weaker k for the weak ones. Now 
we consider a general distribution P(k). We propose the following family 
of variational arguments to define an effective coordination and evaluate the 
density of states. Let us introduce an artificial cut off A > 0, such that the 
coordination number z(X) = z c + A. We use this cut-off to decompose the 
interactions in two groups: the z(X)/2 strongest ones, which form a network 
of " strong" interactions of coordination z(X), and the rest. This allows us 
to repeat the argument of the previous section: we first compute the energy 
u* 2 (X) of the low- frequency anomalous modes of the "strong" network. Then 
we evaluate the correction in energy ((X) induced by the rest of the interac- 
tion on these anomalous modes, as computed in the last section. Finally, one 
obtains the energy E(X) = uj* 2 (X) + £(A) which gives a bound of the energy 
of appearance of the excess-modes. Then the standard procedure consists 
in minimizing E(X) with respect to A. If Aq minimizes this quantity, then 



JE(Xo) gives the best estimate of the frequency of the anomalous modes. 
Furthermore, one can define an effective coordination z = z c + Ao- The 
length that characterizes the anomalous modes follows, according to Chap- 
ter HI I* ~ Aq ■ In what follows we compute E(X) and Ao, discuss the quality 
of the estimate obtained by this argument, and compute other quantities 
that should enable to test if the excess-modes of Lennard- Jones systems are 
correctly described in terms of anomalous modes. 

We first evaluate u;* 2 (A). In the harmonic case where all the stiffnesses 
are 1, we had u* 2 = Ai5z, where A\ can be deduced numerically from the 
data of [50J. From Fig l4.ll one has in three dimensions A\ ~ 0.12. In the 
network of coordination z(X), the amplitude of the stiffnesses fluctuates. To 
present this argument in its simplest form, we neglect these fluctuations, and 
replace the stiffness of each contact by the average stiffness as we discussed 
in the first part of section 18.21 We also showed there how a more accurate 
argument could be made, and we shall come back to this issue later. From 
our approximation follows that SE ~ (k), and in our notation we obtain 
lu* 2 (X) = A 2 (k) x X 2 , where (k) x = ^ J^ x) kP(k)dk. jfe(A) is the weakest 
stiffness of the network of "strong" interactions, defined as: 





(8.15) 



We now compute E(X) and Aq. From Eq. (j8.14j) we get for the correction 
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in energy £: £(A) = h So P(k)kdk, and finally: 

1 r k W 

E(X) = Al(k) x \ 2 + - P{k)kdk (8.16) 

Ct J 

The effective extra-coordination Ao is defined as dE J^ \\=\ = 0. This equa- 
tion defines a non-trivial minimum in general, since (i) dE d ^ | a=o — — k(0) < 
as required by mechanical stability, and (ii) if the potential decays reasonably 
fast with distance, E(X) ~ (k)\X 2 ~ A — > oo as A — > oo. Using Eq. (|8.15|) . 
the minimization of -E"(A) leads to: 

2A^(fc) Ao + ^A^ = lMAo) (8.17) 

A necessary condition for this variational argument to be relevant is that the 
solution Ao of Eq. (|8.17J) is small (at least smaller than the end of validity of 
the scaling of Fig ]4.1[ say Ao < 2.5 in three dimensions). Assuming that it is 
the case (and checking it self-consistently later), one may neglect the term 
quadratic in A 2 in Eq. (j8.17J) to find: 

x 1 MAo) /010 , 

An 



AdAf (k) 



Ao 



In three dimension using the value of A\ one gets Ao ~ . Thus a first 

requirement for this variational argument to apply is that fc(Ao) must be 
several times smaller than (k)\ . A second requirement is that the bound 
of the energy E(X ) must be relatively small in comparison with the Debye 
energy E D = uj 2 d \ if this bound is too large, our variational argument does 
not lead to a correct estimate of the density of states, and the vibrational 
modes are not well described by anomalous modes. This implies that the 
correction in energy ((X) must not be too large. It is clear that these two 
conditions are satisfied when the interaction potential decays fast enough. 
Consider for example a potential of the form V(r) = r~@. If the potential 
decays slowly, (3 > is small, then the two requirements we just evoked 
are not satisfied. In three dimensions for example there are on average 12 
neighbors surrounding a particle. If j3 is small they all interact with a similar 
stiffness with the particle lying at the center, and one finds that k(X ) is close 
to (k)x , and not several times smaller. Furthermore the interactions with 
the second, third, etc... neighbors are not negligible either, and ( is large. 
Thus our variational argument is not relevant in this case, and we expect the 
vibrations of such system to correspond to plane waves, and not to display 
excess-modes. If (3 is large enough, then the two requirements are satisfied if 
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the system is amorphous. To show that we assume that g(r) does not evolve 
much with the potential considered, as it is generally observed with radial 
interactions in the glass phase. Then, when f3 becomes large, the hierarchy 
among the stiffnesses diverges: if one particles has 2 neighbors at distance 
T\ and r 2 with T\ < r 2 , the ratio of the stiffnesses of these two interactions 
is |j = (^■) /3+2 which diverges when 3 increases. It follows that £ becomes 

negligible, and that A ~ yrr^ — > as 3 —> oo. 

To decide whether or not the anomalous modes that appear in this vari- 
ational argument are responsible for the excess-modes observed in Lennard- 
Jones systems, such as those of jH], one may consider a more precise ob- 
servable that characterizes the anomalous modes. The quadratic energy 
of a mode is the sum the energy condensed in every contact 5E(kij) = 
\kij[(5Ri — SRj) ■ rlij} 2 - Consider the average energy condensed in a con- 
tact of stiffness k, that we denote (5E(k)), where the average is taken on 
all the contacts of similar stiffness k. According to our variational argu- 
ment, for the lowest-frequency anomalous modes (5E(k)) varies as follow: if 
k < fc(Ao), the displacement of the particles in contact are not correlated and 
(5E(k)) = jjk. If k > fc(Ao), the interaction belongs to the rigid network. 
Then the relative longitudinal displacement between particles in contact cor- 
responds to the modulation by a sine on a length scale I* ~ Aq 1 , therefore 
(5E(k)) ~ A;Aq. Thus the curve ^ 5E ^ is a step function, which jumps at 
k = k(X ), as represented in Fig. ()8.3|) . In the next paragraph we shall argue 
that this result is not exact and that the step is not sharp, but smooth. In 
any case, if the excess-modes observed in Lennard- Jones systems are related 
to the anomalous modes describe here, they should exhibit this cross-over. 
This could be verified numerically. Furthermore, from this curve k(X ) can 
be characterized: it is the stiffness at which the cross-over takes place. From 
k(X ) one can compute the effective coordination A if the distribution of 
stiffness P(k) is known. 

To conclude, we expect that the present variational argument could be 
improved in several ways. For example the fluctuations of stiffnesses could 
be taken into account to estimate u>*, as discussed at the end of section IHT^l 
This leads to corrections on the quantity (SE(k))/k for k > k(X ) that can be 
estimated using effective medium theory. One finds for the anomalous modes 
that the quantity (8E(k))/k is not constant for k > k(Xo), but decreases with 
k, as represented in the full line curve of Fig. lj8.Hj) . Such effects smooth the 
step function drawn in Fig. (j8.3j) . 
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Fig. 8.3: (5E(k))/ k as defined in the text vs the stiffness k for an anomalous 
mode. The red dashed curve corresponds to the variational argu- 
ment presented in the text. The black continuous curve is what we 
expect qualitatively when the variational argument is improved. 



9. Rigidity of hard sphere liquids near 
the jamming threshold 



There is something mystifying about the properties of a hard sphere system. 
On the one hand it is a simple system: no energy is involved, the only rule 
is that particles cannot interpenetrate, and only entropy matters. On the 
other hand, it displays a rich phenomenology. When the packing fraction <p 
is slowly increased from the liquid phase, it crystallizes. If the <ft is increased 
rapidly, crystallization is avoided, and a glass transition is observed. The 
time r a that characterizes the de-correlations of the density fluctuations at 
some vector q grows rapidly. For a 3-dimensional monodisperse system r a 
becomes inaccessible numerically above 0o ~ 0.59. Nevertheless the packing 
fraction can be increased further until the pressure diverges, when the dis- 
tance between particles vanish: it takes place at <p c ~ 0.64, the random close 
packing. At packing fractions between O and C the structure of the system 
is frozen, apart from the fast rattling of the particles around their average 
position. 

The mode coupling theory furnishes predictions for the relaxation of the 
density fluctuations near the glass transition in rather good agreement with 
empirical data |98j . Nevertheless there is no clear understanding of the spatial 
nature of the events that relax the system, neither of the heterogeneous 
dynamics that has been observed near 0o • A necessary first step to study such 
questions is to understand the cause of the rigidity of a hard sphere system at 
times scales t <^ r a . In the conventional picture jM], the freezing of the liquid 
at times t <C r Q is interpreted with the "cage" effect. As the density increases, 
the cages formed by the neighboring particles tighten, and the characteristic 
time for a particle to escape its cage increases. This description considers the 
stability toward the motion of one single particle. It is dangerous since, as we 
discussed in the previous Chapters, the stability against collective motions 
of particles is more demanding than the stability against individual particles 
displacements. For example in d dimensions d + 1 particles are enough to 
pin one particle. Nevertheless, a system with a coordination number d + 1 is 
unstable to cooperative motion, as shown by the Maxwell criterion. 

In this Chapter we study the microscopic cause of the glass rigidity. In 
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particular we study the hard sphere glass at packing fraction near <p c , and 
we derive some elastic properties of the glass phase. Following the ideas 
of the previous chapters, we show that the solid-like behavior at i <C T" a 
requires the formation of a rigid structure with a sufficiently large coordina- 
tion number. We argue that such structure corresponds to the network that 
carries momentum between the particles, which was introduced recently to 
study granular flows |lUUj . Our main achievement is to compute an effec- 
tive potential between particles in contact through this network. We show 
theoretically, after averaging over the fast fluctuations or rattling of the par- 
ticles, that the hard sphere potential becomes logarithmic. This result is 
exact at C , and is in very good agreement with a direct numerical check. 
This allows us to define normal modes and to apply the results we found for 
soft sphere solids near the jamming threshold. In particular, the extended 
Maxwell criterion applies, which yields an inequality for the network coordi- 
nation. We confirm this result numerically. We compute the scaling of the 
high-frequency elastic moduli near <fr c . More generally, this approach shows 
that the jamming threshold acts as a critical point both in the solid and in 
the glassy liquid phase. This suggests original relaxation processes in the 
super-cooled phase. 



9.1 Coordination number and force 

For concreteness we consider a hard sphere system without dissipation where 
particles collide elastically. We show in the next section how to generalize 
our finding to the dissipative case where particles follow Brownian motion. 
We consider packing fractions such that the typical collision time between 
two neighbors r c is much smaller than r a . We introduce an arbitrary time 
t\, much larger than the collision time, and much shorter than any time 
scales at which the structural relaxation occurs. Two particles are said to 
be in contact at a time t if they collide at least once in the time interval 
[t — ti/2,t + ti/2]. This enables us to define the coordination number z as 
z = 2N C /N, where N c is the total number of contacts. Then, we define the 
contact force as the total momentum per unit time exchanged between 
the two particles: 

i n=n col 

h = T E AP n (9.1) 

'1 n=l 

where the sum is made on the total number of collisions n co \ between % and j 
that took place in the time interval [t— ii/2, t+ii/2], and AP n is the momen- 
tum exchanged at the nth shock. These definitions were first introduced in 
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a work on dense granular flows |lUUj , and were used recently on hard sphere 
systems jZHl HOlj . Coordination and contact forces depend a priori on an 
arbitrary parameter t\. In the high-packing fraction system we studied, we 
did not observe any relevant dependence of these objects with t± as long as 

In Fig (j9.1j) we show a two-dimensional example of the contact force field 
obtained with such procedure at packing fraction <fr very close to <p c . Note 
that the forces are roughly balanced on every particle, as it must be the case 
on time scales over which the structure is stable. 

To obtain high packing fraction configurations such as the one of Fig. 1)9.1)) , 
we proceed as follows: we consider the two-dimensional polydisperse 23 con- 
figurations of [201 at the jamming threshold at packing fraction <p c m 0.83. 
At this packing fraction the particles are in contact. Then we reduce all the 
particles diameter of a relative amount e. This leads to configuration of pack- 
ing fraction <ft = <f) c (l — e) 2 . We assign a random velocity to every particle. 
A Newtownian dynamics is then computed using an event-driven simula- 
tion. As we shall discuss later, such a protocol does not lead to a system 
at thermal equilibrium. Nevertheless, we are not interested in thermody- 
namic properties, and in practice systems with such high packing fraction 
are never equilibrated in any reasonable time. We rather aim to study the 
conditions that guarantee mechanical stability. Such condition should be ful- 
filled whether the system is at thermal equilibrium or not, as long as it is 
stable on reasonably long time scales. 

Note that since there is no energy involved in such system, the temper- 
ature only fixes the time unit. In what follows we impose for the average 
square velocity (v 2 ) = 1. 

9.2 Effective potential 

In the previous Chapters we studied the rigidity of amorphous solids by 
considering their vibrational modes. It is a priori problematic to use the 
same analysis to study the rigidity of contact networks such as the one of 
Fig.(jnHJ): the hard sphere potential is discontinuous, and the energy cannot 
be expanded as in Eq. ()2.7j) . Nevertheless we shall go around this difficulty 
by deriving a smooth effective potential. The trick is to average on the fast 
fluctuations that take place on characteristic time smaller than some t% 3> r c . 

22 If t\ is too small some contacts can disappear, which can lead to the appearance of 
unstable modes when they are computed following the procedure bellow, even thought the 
system is stable 

23 Half of the particles have a diameter unity, the other half have a diameter 1.4. 
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Fig. 9.1: Contact forces for iV = 256, e = 1CT 5 and t\ = 10 5 . Black points 
represent particles. Contact forces are sketch by arrows which 
start from the particle center, and whose length is proportional to 
the force amplitude. Note that the forces are balanced on every 
particle, as it must be the case on time scales where the structure 
is stable. For similar force networks see |101j . 
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Consider two particles i and j in contact separated by a spacing — — rj 
of typical value h, where r\ and Tj are the radii of the particles, and 
the distance between % and j. Such a spacing fluctuates in time between 
(when the particles are colliding) and a few h, so that the instantaneous 
value of this spacing does not give much information about the contact ij. 
Nevertheless, if r\j — r j — r,j is averaged on large time intervals t\ r c , the 
spacing hij = (r\j — r^ — rj) converges to a well defined value. We shall define 

— * 

the averaged position R av {t) of particle i as: 

R7(t) = T f t+th m')dt'. (9.2) 
ti Jt-ti 

In what follows we estimate the contact force exchanged between two 
particles % and j whose spacing is ~ \\Rf v — Rj V \ \ — n — rj. We furnish 
thermodynamic arguments, that apply both to Newtonian and Brownian 
dynamics. We start by considering a one-dimensional system with beads of 
diameter 1. The partition function Z is: 

nrhi=oo 
/ dh ie - ph > (9.3) 

i Jhi=0 

where hi is the spacing between particle i and i + If an external force dipole 
Pi = —pi+i = Pi is applied on i and i + the partition function becomes: 

Z = \\ dh je - ph i / dh,e~ {p+p ^ (9.4) 

^A 3 =0 Jhi=Q 

From the partition function one can compute the average spacing (hi) = ^f^-- 
Since the contact force /j in the contact i, i + 1 is — p± + p, one finds: 

l = \ (9.5) 

This result can be extended to an isostatic state of any spatial dimen- 
sion. A particularity of the isostatic state is that the number of displacements 
degrees of freedom is precisely equal to the number of contact. Hence the 
configuration of the system can be described by the set of distances between 
particles in contact. If the system is at equilibrium in a meta-stable state 
where the contact forces field |f) = {fij} is well-defined, the partition func- 
tion can be written: 

rhij=oc 

Z = \{ dhijcr^'"! (9.6) 

(ij) Jh i3= G 
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The argument valid for the one- dimensional line of particles is valid here, 
and one obtains (hij) = f^ 1 . This shows that in an isostatic assembly of 
colliding hard spheres, when the particles rattling are averaged, the hard 
sphere potential converge to an effective potential. At time t>r c the system 
behaves as an assembly of particles of positions \H av ) = {-R™} interacting 
with the potential V^r): 

Vij(r) = 00 if r < r j + r j 

Vij(r) ~ —ln{r — — rj) if % and j are in "contact" 

^'( r ) = if i and j are not in "contact" (9.7) 

The relation force/ distance is checked numerically in Fig. ()9.2j) at packing 
fraction close to C . At such packing fractions the system is nearly isostatic, 
as we shall see in the next section. The exponent found is in very good 
agreement with Eq. (J9.5|) . Fig. ()9.2|) also shows that for large ti the disper- 
sion of the contact forces around their average value described by Eq. (j9.5|) 
is extremely small. This indicates that the only relevant parameter that 
characterizes the contact force amplitude is the spacing h, as predicted by 
Eq.flHH). 
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h 

Fig. 9.2: Log- log plot of the contact force amplitude versus the spacing h = 
r — Ti — rj for various e in a system of N = 256 particles. Each dot 
represents the pair of number (fij, {hij}) associated with the contact 
ij. Dots collapse on the dotted theoretical curve defined by Eq. (|9.5j) . 

When the coordination z increases, the hij are not independent variables 
anymore. We aim to evaluate the corrections to Eq. (J9.7|) when the system 
is at a finite distance from isostaticity, that is with 5z > 0. Isolating the 
contact ij we may write the partition function as follows: 

Z= / dhtje-te^e-*^ (9.8) 

J hij =0 
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where T(hij) is the free energy of the entire system conditionally to the 
value of the spacing h%j. It can be written as Tihij) = J-'o + A J 7 (hij), where 
J-'o does not depend on hij. To evaluate AJ-"(hij), we consider as a zero 
order approximation that Eq. ()9.7j) is true. Then AT(hij) corresponds to the 
energy cost induced by a local strain of amplitude — (hy) = hy — f^ 1 
in an elastic system where particles interact with the potential (|9.7jl . In 
Chapter [7| we computed the response to a local strain of amplitude e. The 
corresponding energy cost varies with the contact considered, and for small 
strain its amplitude follows SE ~ 5zBe 2 , where B is the bulk modulus. For 
an interaction given by Eq. ()9.7|) . following chapter [7| one finds B ~ (V"(r)) ~ 
((r — 7-j — rj)~ 2 ) ~ p 2 , so that AJ-(hij) ~ 5zp 2 (hij — fij 1 ) 2 - Thus we may 
write A.T{hij) = bzC^f 2 j{h^ — fij 1 ) 2 , where CV,- is positive, of order one, 
and can a priori depend on the contact considered. Using this expression in 
Eq. (|9.8|) . and expanding Z to first order in 5z, we can compute the corrections 
to (hij). One finds (hij) = -r-(l — 2C i j5z) ) so that the force-displacement 
relation satisfies: 

f ij = 1 ^-(l-2C i j5z) (9.9) 
hij 

This estimates the corrections to the potential of Eq. (j9.7|) . which vanish 
when the system becomes isostatic. Thus, one non trivial consequence of 
these corrections is to weaken the force for a given inter-particle distance. 
To test this effect we compute numerically C(Sz) = (fij (hij)) ij — 1, where 
Qij denotes the average over all contacts. The results are represented in 
Fig. ()9.3|) . Small corrections are indeed found, which are in good agreement 
with our prediction C(5z) ~ — Sz. In what follows we are mainly interested 
in scaling relations near C , where corrections of the order 5z do not matter. 
Thus we shall neglect them, and consider that the effective interaction is 
constant and given by Eq. (j9.7|) . 



Note that Eq. flEfl 

can be recovered with a simple scaling argument. We 
may evaluate the collision frequency v as v ~ where v is the average 
amplitude velocity. The typical momentum exchanges during one collision 
is, taking m = 1 for the particles mass, ||AP|| ~ v. Therefore following 
Eq. (j9.1|) we recover that the force follows /y ~ ~ 

In what follows we shall only assume that the potential V(r) is differ- 
entiate and that the force f(h) scales as h^ 1 . The corrections estimated 
above in Eq. (j9.9|) do not affect these results, therefore we neglect them and 
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Fig. 9.3: Average correction C{5z) as denned in the text vs. excess coordination 
Sz for various (p. The line is a linear fit consistent with our predictions 
at small Sz. Corrections are small, of the order of 3-4 percents when 
Sz = 1. 



use Eq. ([9.7|) . As we shall see, this implies that <ft c acts also as a critical point 
in the liquid phase. In particular, Eq. (j9.7j) enables us to define the stiffness 
kij of a contact ij as: 

k " = v » ~ v^hw ~ % (910) 

This allows to define a dynamical matrix and vibrational modes once the 
average particles positions and the contacts are known. As we discussed 
above, the potential of Eq. (j9.7j) has an entropic nature. Thus such vibrational 
modes describe the local curvatures of the entropic landscape of the system. 
In what follows we show that imposing the stability of these modes leads to 
constraint on the coordination of the force network. Then we discuss the 
elastic property of hard sphere glass, derive the elastic moduli and discuss 
the length scales that appear in the response of such systems. Finally we 
discuss how these modes may be related to the structural relaxation. 



9.3 Stability of hard sphere systems 

If the network of contact of a hard sphere system is weakly connected, we can 
apply the results of Chapters El Eland El In particular, following Eq. ([5.5p and 
Eq. (j8.11|) and the paragraph above it, we obtain that such system presents 
anomalous modes that appear at a frequency oj 2 am = AiB{p)5z 2 — A 2 p, 
where B(p) is the bulk modulus of the system. The bulk modulus scales 
as the average stiffness of the contacts, and following Eq. (|9.1()jl we obtain 
B(p) ~ p 2 . As we discussed in Chapter a rigid system does not display 
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Fig. 9.4: Density of vibrational modes D(u) versus frequency in a hexagonal 
crystal of 1024 particles for e = 10~ 4 . The frequencies are rescaled 
by e _1 . The particle positions were average in time t\ — 4 x 10 4 to 
obtain |R av ). Eq. ()9.7|) was used to compute the dynamical matrix, 

from which the vibrational frequencies were inferred. 

** 

unstable modes. Therefore we obtain that there must be a constant Co such 
that: 

5z > C p'^ (9.11) 

which is another realization of the extended Maxwell criterion we derived in 
Chapter Note that near the jamming threshold, the typical inter-particle 
spacing e goes as e ~ C — and thus p ~ e _1 ~ (0 C — as was already 
computed with different methods jZHlEO]- Thus Eq. (j9.11|) indicates that 5z 
must scale as Sz ~ (<p c — cf)) 13 with (3 < ^, as is the case for a soft sphere 
system above 4> c . 

9.3.1 Stability of the hexagonal and the square 
crystals 

To test Eq. (j9.11jl we perform tests on three different systems in two dimen- 
sions. We start by considering an hexagonal monodisperse crystal and a 
monodisperse square crystal. We consider these two systems at their maxi- 
mum packing fraction where hard spheres are in contact. Then we reduce the 
particles diameter by a relative amount e, and start the dynamics. According 
to Eq. (j9.1H) these systems must behave very differently. In the hexagonal 
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crystal, the coordination is 6, therefore 5z = 2 3> ~ ea for small e. 
Therefore the condition of Eq. (j9.11j) is satisfied and we expect the system 
to be stable. On the other hand, in the square crystal, the number of first 
neighbors is 4, 5z = 0, and the system cannot satisfy Eq. (j9.11j) without large 
structural rearrangements for any e. Thus such a system cannot be rigid. 
These predictions are verified numerically. For small e, the hexagonal crystal 
displays no structural changes, whereas the square crystal collapses rapidly. 
Such collapse leads to an hexagonal configuration, and is generated by the 
buckling of unstable modes. We show the corresponding displacements in 
Fig. fj9.5)) . The instability can be also observed in the free energy landscape 
by computing the vibrational modes of the corresponding systems. Fig. (j9.4j) 
shows the density of vibrational modes of the hexagonal crystal, which varies 
linearly at low frequency, as expected for a two-dimensional crystal. No un- 
stable mode are observed. In the square crystal, see Fig. (j9.(ij) . the density 
does not vanish at u — > 0, as expected for an isostatic system. Furthermore, 
we observe unstable modes, as implied by Eq. (j9.11|) . 




9.3.2 Stability of hard sphere systems near the 
jamming threshold 

We perform the same test using the polydisperse configurations of [SU] at 
the jamming threshold. We study the structural stability of the system by 
computing the self correlation function C(q,t) = (exp[iq ■ (Rj(t) — Rj(0))]) 
where the average is taken on every particle j. In what follows we consider 
wave vectors q of norm where r\ is the radius of the smallest particle. 
Typical curves for different e are represented in Fig. (j9.7|) . For e smaller that 
roughly e = 0.05, the system ages and dynamics is intermittent. There are 
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Fig. 9.6: Same caption as Fig. ([9.4)1 for a N=1024 square crystal. 80 
unstable modes were observed, lying in the frequencies range 
uj G [—0.015,0]. There are very degenerated and therefore not 
appropriate for a plot, thus we do not represent them. 




Fig. 9.7: Examples of C(t) = C(q = l,t) for q = -njr\ versus real time for 
different e in a N = 256 particles system, t = corresponds to the 
initial time where random velocities are assigned to the particles. 
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long periods of time where the structure is stable, that appear as plateaus 
in the self correlation function. Such periods ends with crashes, where the 
self correlation function drops of a large amount in a very short time. These 
events must correspond to sudden collective rearrangements involving a large 
number of particles. During such crashes, we observe that the coordination 
increases, and that the pressure drops. We also observed that the correlation 
function of the force network H(t) = (fij{T)fij(r + t)) is constant during the 
plateau, and drops when a crash occurs. 

Such crashes, or "earthquakes", were reported in other glassy systems. 
They occur when a Lennard- Jones liquid is rapidly quenched at temperature 
much lower than the glass transition |lU2j . In this case the crashes typically 
involve 100 particles. Such crashes were also observed in colloidal pastes using 
dynamic light scattering |10Hj . and in the dielectric response of laponite [104] . 
It would be obviously interesting to understand what is the nature of these 
crashes, and what triggers the sudden collapse of apparently mechanically 
stable structure. These are fundamental questions of the glass transition. 
In the present Chapter we study the stable structures that appear before 
and after the crashes, when the system is quiet. Understanding why such 
structures are rigid is certainly a necessary step to find out why and how 
they can yield. 

We define the quiet periods as the plateaus of C(t). To study the rigidity 
of the dense hard sphere assemblies during these periods, we compute the 
vibrational modes. If the rattlers are removed 24 , we find that, the system 
is mechanically stable: there are no unstable modes. The density of states 
for e = 10 -4 is shown in Fig. ()9.8j) . The main difference with the other stable 
structure that we studied above — the hexagonal crystal — is the presence 
of a large excess of modes at low frequency. A large amount of modes are 
nearly unstable, as we expect near the jamming threshold where we must have 
D(lo) — > u° when the frequency vanishes. Note that the density of states does 
not have a flat plateau as the harmonic soft spheres at the isostatic point, 
despite that the coordination network is the same when e is small. The main 
difference between these systems is the stiffness disparity: in the amorphous 
hard sphere system, the disparity of the contact force leads to a disparity in 
the stiffness since k ~ f 2 . If this disparity is removed by keeping the same 
contact network, but by imposing a constant stiffness in all contacts, we find 
that a flat plateau is recovered. 

The absence of unstable modes enables us to test Eq. (j9.11j) . that must be 

24 The rattlers do not participate to the rigidity of the structure. They can be identified 
as the distance with their neighbors is much larger than for the rest of the inter-particle 
distances, so that the frequency of the shocks they have with their neighbors is much 
smaller than the other particles. 
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Fig. 9.8: Density of vibrational modes D(u) versus frequency in a amor- 
phous hard sphere glass of 256 particles for e = 10~ 4 . The fre- 
quencies are rescaled by e _1 . D(uS) was computed in a quiet period 
preceding the first crash. No unstable modes were observed. 

satisfied as the system is stable. To check this relation we computed numer- 
ically both the coordination and the pressure for various packing fractions, 
and for various stable periods that appear along the aging regime. As shown 
in Fig. (jHZIl) , the data are consistent with an equality of the inequality 1)9.1 lj) . 
This suggests that a hard sphere glass is only marginally stable, as it is the 
case for soft spheres slowly decompressed toward <p c . 

9.4 Elastic property of the hard sphere glass 

We now use this approach to derive the elastic behavior of the hard sphere 
glass near C . The results of Chapter apply. In this system the elastic 
moduli have purely entropic causes: they describe how the number of con- 
figurations is reduced under a global strain. As we already discussed, the 
scaling of the bulk modulus B is: 

B ~ p 2 ~ (4> c - (f))- 2 (9.12) 

In repulsive systems near C the shear modulus scales as the bulk modu- 
lus times the excess coordination. Thus the extension of Eq. ()7.HHJ) to non 
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Fig. 9.9: Log-log plot of 5z versus the average contact force (/) ~ p. The data 
were obtained for different e and different time periods. The black line 
corresponds to the equality of the inequality (|9.11j) . 



harmonic contacts yields: 



G ~ B5z ~ 5zp 2 



(9.13) 



If one assumes that the glass is marginally stable, that is to say that the 
system lies on the bound of inequality (|9.11j) . one obtains: 



G 



P 



3/2 



(9.14) 



As in practice such glass is not at equilibrium near <p c , the hypothesis of 
marginal stability certainly depends on the preparation of the system. It 
is reasonable to think, following the analogy with the soft spheres, that if 
the system is slowly compacted toward C the dynamics would not bring the 
system far away from marginal stability. Thus we expect 8z ~ to be a 
good approximation and the shear modulus G not to differ to much from the 
scaling law of Eq. ()9.14|) . 

Furthermore, this system present excess modes: there are the anomalous 
modes we discussed in the previous Chapters. Near the jamming threshold, 
when the pressure diverges, we expect D{oj) — > u°. When cf> decreases, the 
density of states is characterized by some frequency uj* ~ B^{p)5z ~ pr^bz. 
Consequently a hard sphere glass is characterized the length I* ~ 5z~~ l which 
characterizes the anomalous modes present at low-frequency. 



9.5 Discussion 



In this Chapter we studied some properties of the hard-sphere glass near 
the close packing at <p c . We argued that the mechanical stability at time 
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t <C r a relies on the emergence of a rigid force network. The key point 
was to compute the properties of the force networks by averaging on fast 
fluctuations, the particles rattling. This allows one to derive an effective 
potential, exact at the isostatic limit, and to define rigorously the vibrational 
modes in hard sphere systems. They characterize the entropic landscape 
expansion around a configuration. This allows one to recover the main results 
valid for soft sphere systems near the jamming threshold. In particular the 
stability of the vibrational modes imposes a bound on the coordination of the 
contact force network. More generally, we showed that <p c acts as a critical 
point: the elastic moduli and the coordination z scale as — > C . The system 
responses are characterized by I* that diverges at C . 

As critical points and diverging length scales are in general associated 
with diverging time scales, the presence of a critical behavior at <j) c may be 
the cause for the glass transition observed in hard sphere systems. In partic- 
ular, our work suggests that the slowing-down of the dynamics is due to the 
appearance of an extended rigid network of interactions, rather than to any 
local properties such as the tightening of cages. This suggests an alternative 
perspective for the relaxation in super-cooled liquid, at equilibrium or in an 
aging regime such as the one of Fig. (j9.7j) . The particularity of the systems 
that crash, such as the square crystal or the amorphous state near isostatic- 
ity, is the large amount of modes around zero frequency. As we discussed in 
previous Chapters such modes are very different from plane waves, as they 
have rapid spatial fluctuations, which makes them sensitive to the applied 
stress term responsible for buckling. Furthermore they can be localized on 
length scales larger than some /* without changing in their frequency. These 
modes indicate a nearby elastic instability, and may well play a role in the 
rigidity loss. Furthermore, whatever rearrangements occur when the system 
relaxes, suddenly or not, the corresponding displacements must be similar to 
the soft modes, as there are, rougly speaking, the only modes where particles 
can avoid to inter-penetrate. Thus it is tempting to associate the structural 
relaxation with the structural buckling of the weakest frequency modes, that 
could be induced for example by pressure or coordination fluctuations. If so, 
we expect the relaxation to occur on length scales of order I* ~ pz, since 
for smaller subsystems the frequency of the modes increases, and the elas- 
tic instability goes away. Thus our work suggests an original perspective on 
the possible cause of the heterogeneous dynamics that occur in super-cooled 
liquids [2U 120] • Many models of the glass transition lead to heterogeneous 
dynamics [221 ■ in most of them, a local rule describes the motion of one 
or few particles, and leads to cooperative dynamics. Our work suggests an 
alternative view: the rule that allows the particle motion is itself non lo- 
cal, for the simple reason that rigidity is not global property, as was already 
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understood by Maxwell more than one century ago. 



10. Conclusion 



10.1 Summary 

The elastic properties of an assembly of short-range, repulsive particles dis- 
play critical behaviors near the jamming transition. In particular the elastic 
moduli, the frequency of the excess-modes, and the coordination number 
scale. In the first part of this thesis we derived the corresponding exponents. 
Then we showed how some of these results can apply to various systems, 
such as silica glass or colloidal particles. 

The starting point is the following variational argument: if a rigid system 
of coordination z is cut in sufficiently small subsystems of size /, these sub- 
systems are not rigid. This is true as long as the subsystems size is smaller 
than some length I* ~ 5z~ x that diverges at the jamming threshold. Smaller 
subsystems contain modes of zero frequency, the soft modes. From these 
modes one can build what we called the "anomalous modes" which have a 
frequency of order 1// in the original system. This gives for the dependence 
of the onset of excess-modes oo* ~ 5z, as observed numerically. The system 
can be described as a continuous elastic medium only for uj < uj*. At larger 
frequencies, it behaves as an isostatic state. 

Then, by arguing that the soft modes have large transverse relative dis- 
placements, we showed that the anomalous modes were much more sensitive 
to the applied stress than the conventional acoustic modes. This has two 
direct consequences: in a purely repulsive system, anomalous modes can ap- 
pear at frequencies much smaller than uj* . Furthermore, in rigid solids the 
coordination number must be sufficiently large to guarantee the stability of 
the anomalous modes. We find Sz > Co{p/ B)^ ~ (0 — C ) 5 , where p is the 
pressure and B the bulk modulus. This generalizes the Maxwell criterion 
for rigidity 5z > valid when applied stress is absent. It follows that, if 
the jamming threshold C is reached adiabatically, the pair correlation g(r) 
measured at C must display a divergence g(r) ~ ( r -i)i w ^ n 7 — \- This 
divergence is the vestige of the excess contacts 5z necessary to maintain the 
rigidity of the structure at larger 0. 

These arguments furnish the scalings of both the onset of excess modes 
and the coordination number, but they do not enable us to compute the elas- 
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tic moduli. To do so, we used the fact that the linear equation that defines 
the soft modes is the dual of the force balance equation. This allows us to 
derive an original formalism relating the response to a strain to the contact 
force fields that satisfy force balance on each particle. This enables us to 
compute the shear and the bulk modulus. We found that for a repulsive 
system the ratio G/B vanishes at the jamming transition. We also obtain 
exact results on the response to a local strain. The energy cost of such de- 
formation vanishes at the transition as the shear modulus. At the transition, 
this response extends in the entire system. Therefore this is also true for 
the soft modes that appear when one contact is cut in an isostatic state, an 
assumption that was essential to compute the frequency of the anomalous 
modes. 

In a second part we studied the applications of these concepts to glasses, 
granular matter and colloids. All glasses have an excess of low-frequency 
modes, the boson peak. It is especially strong in silica, one of the best 
glass former. In silica the strongest interactions are those responsible for the 
rigidity of the Si0 4 tetrahedra. If the weaker interactions, such as Van de 
Waals, are neglected, the tetrahedral network obtained is isostatic. Following 
our argument, such a network must have a constant density of states at low 
frequency, as it is observed numerically. When the weaker interactions are 
turned on, the anomalous modes responsible for the plateau shift to higher 
frequencies. This improved variational argument predict the appearance of 
a plateau in the density of states. Such plateau is indeed observed in the 
simulations of silica, and appears around 1 THz, the boson peak frequency. 
As the key parameter of our description is coordination rather than disorder, 
this argument also justifies the similarity between the density of states of 
silica and the one of the corresponding crystals, which also present a plateau 
at roughly the same frequency. This is certainly a strong point, as most 
of the existing boson peak theories are based on disorder only, and cannot 
explain the excess modes that show up in these crystals. Finally we proposed 
to extend these ideas to Lennard- Jones systems. 

Generally critical points display similar behavior on both sides of the 
transition. We argued that it is also the case for the jamming transition. We 
showed that when a hard sphere liquid approaches the jamming transition, 
the contact force network, that now characterizes how particles exchange 
momentum, is very similar to the one of elastic spheres above C . The key 
point was to show that when the dynamics of an isostatic system is averaged 
on short time scales, the interactions among particles can be described by a 
logarithmic effective potential. We evaluated the corrections of this potential 
when the coordination increases. This allows to compute the normal modes 
of such systems. As a consequence the results of elastic spheres solids apply 
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to hard sphere liquids. In particular the extended Maxwell criterion must 
be satisfied near <ft c , when the system is rigid on short time scales. This ap- 
proach also yields the elastic moduli that characterize the system for t <C T a , 
the relaxation time. This description supports that the relaxation in super- 
cooled liquid should not be described in terms of the motion of individual 
particles, even if models with such local rules can lead to non-trivial cooper- 
ative dynamics. Rather, it suggests that the elementary motions to consider 
are themselves collective, and are related to the anomalous modes introduced 
above. 

10.2 Perspectives 

10.2.1 Low temperature glass properties 

The transport properties of glasses, but also of granular matter [1051 HU6j , 
are not understood In particular there is a plateau in the thermal con- 
ductivity around 10 K, temperatures at which the heat transport can be 
dramatically smaller than in the crystal phase j^]. These temperatures cor- 
respond to frequencies in the Thz range, where the boson peak appears. In 
repulsive, short-range systems we showed that the cause for the boson peak 
is the following: above some frequency u* the system behaves as an isostatic 
state, whose density of states is much larger than the one of the crystal at low 
frequency. We argued that these ideas are more general and apply for exam- 
ple to silica and granular matter. This suggests that the transport properties 
of glasses in the Thz range correspond to those of a system at the jamming 
threshold. Then the question is to understand how the anomalous modes we 
introduced contribute to the transport. Accordingly, it is necessary to com- 
pute their spatial power spectrum E(q). For a normal mode of displacement 
{SRi}, E(q) is defined as: 



For an acoustic mode of frequency uj, E(q) presents a peak for q = c _1 cj, 
where c is the longitudinal velocity of sound. If this peak has a finite width 
Aq, the scattering length of the acoustic mode follows I ~ Ag _1 . This length 
enters in the computation of the thermal conductivity: the contribution of a 
mode to the heat transport goes as l-c PQ. The empirical data suggests that at 
frequency smaller than the boson peak, the acoustic modes have a scattering 




(10.1) 



(10.2) 
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length I m 150A, where A is the wavelength of the corresponding mode [3]. 
We expect that the anomalous modes that appear at higher frequencies will 
transport much less than that, because their spatial correlations (SR(0) ■ 
5R(x)), which enter in Eq. QlU.lJl . are presumably very small. This quantity 
is similar to the spatial correlations of the soft modes. If these correlations 
were zero, both the soft modes and the anomalous modes would be equally 
distributed on all wave vectors. Then the anomalous modes would have a 
purely diffusive behavior, corresponding to a scattering length / = 1 particle 
size. Thus they would almost not contribute to the transport at all. It is 
apparent from Fig. (J3. 1|) that the spatial correlations of the soft modes are 
indeed small. This observation may furnish an explanation for the observed 
bad quality of the transport at the boson peak frequency. It is also consistent 
with the empirical observation that the boson peak frequency depends only 
weakly on wave vectors q if at all |lU7j . indicating that the excess-modes have 
a wide distribution over the wave vectors q. For a quantitative discussion it 
would be of great interest to derive the spatial correlations induced by the 
soft mode equation 1)2. Furthermore, other interesting effects could in 
principle affect the transport. In particular, our derivation of the anomalous 
modes does not preclude the presence of underlying acoustic modes, that 
could hybridize with the anomalous modes, and enhance the transport. 

At lower temperature, the properties of glasses are still a challenge to the- 
ory. The specific heat has a nearly linear dependence with temperature, and 
the thermal conductivity varies quadratically. This is in general interpreted 
by the presence of two-levels systems: atoms or groups of atoms can switch 
between two configurations by tunnel effect. This model is phenomenolog- 
ical and there is no consensus on what these two level systems may be. It 
has been often argued that the excess-modes of the boson peak are good 
candidates to form two-levels systems, see e.g. jUSUHHl- As these modes are 
soft, non-harmonic terms are important, which could lead to the canonical 
form of 2-levels systems: two wells separated by a potential barrier. If it 
is so, our interpretation of the boson peak suggests that the 2-level systems 
are rather extended, plausibly on the length I*, and that the displacement of 
each particle could be much smaller than a particle size. It should be possible 
in principle to test this possibility. As we discussed earlier, the pressure has 
two opposite effects on the anomalous modes: on the one hand it increases 
the coordination, and on the other hand the applied stress term lowers their 
frequency. We expect that in some systems, the destabilizing effect of pres- 
sure dominates. This can even lead to an elastic instability where anomalous 
modes become unstable, and where the configuration of the system changes. 
Such elastic instability might occur for example in silica glass at high pressure 
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[108 . If the density of anomalous modes is increased at very low temperature 
by tuning the pressure, the density of two-levels systems should be enhanced 
too. The latter could be directly checked by specific heat measurements. 

10.2.2 The glass transition 

We discuss the following possibility: in fragile glass former a fast slow down 
of the dynamics occurs close to the temperature at which the system manages 
to form an extended rigid network. We first study the role of temperature 
on the anomalous modes in the glass phase. Then we propose a microscopic 
distinction between strong and fragile glasses 26 , and we discuss a possible 
structural relaxation process in the fragile case. 

Role of temperature 

Following the results of Chapter |H1 we can write for the onset of the anoma- 
lous modes at zero temperature ui 2 AM = A\Bbz 2 — TC(p) + C, where B is the 
bulk modulus, 5z is the effective excess-coordination number, ( quantifies the 
effect of the weak interactions on the anomalous modes, and Ai is a numer- 
ical constant. 7i(p) is the correction induced by initial stress term. For soft 
spheres near the jamming threshold where the distances between particles 
in contact are similar we found 7i ~ A\p ~ (/), where (/) is the average 
contact force. In a system where the distance r between interacting particles 
can vary, following Eq. fj2.9j) we have Ti, ~ (/ /r). In the present qualitative 
discussion we shall neglect these corrections and consider Ti ~ Aip. 

If the system is heated at constant pressure, we may extend this equation 
and write: 

u 2 AM (T) = A.BiT^ziT) 2 - A 2 p + C(T) (10.3) 

In most glasses, when the temperature increases at constant pressure, the vol- 
ume grows, and B decreases. As the typical inter-particle distances grows, 
we expect both the effective coordination and the effect of the weak interac- 
tions to diminish. Thus all the positive terms in Eq. ljlO.Hj) decrease. On the 
other hand, the pressure is constant. Hence u AM (T) decreases, as is indeed 
observed in most glasses. Eventually u AM (T) reaches zero frequency, as it 
has been observed numerically |J2] and empirically j2H] ■ We shall denote the 
temperature at which rigidity is lost T r . Note that in few materials, such as 
silica, the volume decreases and the bulk modulus grows with temperature. 

26 The relaxation time of strong glass as an Arrhenius dependence with temperature, 
and grows faster in a fragile glass. 
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In this case Eq. fjl0.3j) predicts that the boson peak shifts to higher frequency, 
as observed empirically, see e.g. [IE]. 

Note that when T > T r , there are continuously unstable modes. Interest- 
ingly, such a temperature also exists in mean field spin models 1 1 00' . which 
were proposed as possible scenarios of the glass transition |110j . The dy- 
namics of these models at higher temperature is exactly described by mode 
coupling equations According to this analogy T r = T MC t- 

Fragile and strong glasses 

In order to discuss the distinction between strong and fragile glasses, we 
introduce a second characteristic temperature T a . T a corresponds to the 
typical energy activation of the structural relaxation processes that do not 
depend on the rigidity of the structure. In particular we think about local 
rearrangements, such as the displacements of coordination defects in strongly 
covalent networks [112J. The glass transition must occur around the smallest 
of the two temperature T a and T r , since by definition the structure relaxes 
easily at higher temperature. Hence if T a 3> T r the glass transition occurs 
in the vicinity of T r . As the curvatures of the energy landscape dramatically 
evolves with temperature near T r , it is natural to expect the corresponding 
dynamics to be super-activated. This supports the following scenario: glasses 
with T a ^> T r are fragile. On the other hand, if T a <C T r , the glass transition 
takes place in the vicinity of T a . As the local structure does not display 
any important changes near the glass transition, we expect such glasses to 
be strong. A similar discussion in terms of energy landscape is presented in 

This scenario implies that strong glasses are whether (i) system with a 
large coordination, where anomalous modes frequencies are high. This is 
coherent with the empirical fact that strongly connected covalent networks 
are strong, (ii) anomalous systems which have a strong boson peak, but where 
an increase of temperature does not lower much the boson peak frequency, 
or even stabilizes the system like in silica. In these two cases we expect 
the rigidity of the covalent network not to do play any role at the glass 
transition. This is supported by simulations that indicate that the structural 
relaxation is purely local in silica glass |112j . and that the covalent network 
exists until 8000 K at our pressure |113j . which is much larger than the 
glass transition temperature. Note the elastic anomaly of silica disappears 
at high pressure, which suggests that this glass might become fragile when 
the pressure increases |115j . 

According to the present point of view, fragile glasses must display anoma- 
lous modes that approach zero frequency around the glass transition. A pos- 
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sible test could be done by considering particles with gaussian potentials. It 
was shown theoretically and numerically [J2j 0H1 HI] that such systems at in- 
finite temperature display excess-modes. These modes become stable above 
a finite density. Above this density this scenario predicts that the glass is 
strong, as the stable excess-modes at infinite temperature will stay stable as 
the temperature is lowered. At smaller density we expect a fragile behavior 
to occur, as it is the case for soft spheres |117j . 

Heterogeneous relaxation in fragile glass 

According to the present scenario, in a fragile glass at T r the anomalous 
modes are characterized by a finite length scale I*. According to Eq. fjlO.Mj ) 
/* ~ bz~ x does not diverge since the initial stress term, and the pressure, have 
a finite value. Thus we also expect the shear modulus at T T to have a finite 
value, as it is the case for the marginally rigid system of jSHj- As soon as 
T >T r , the system displays unstable normal modes. Near T r , only the modes 
with characteristic length I* are unstable. The anomalous modes confined 
on subsystems smaller than I* have higher, non-zero frequencies. Thus the 
collapse of unstable modes at T r involve rearrangements on length scale of the 
order of /*, but not smaller. As the temperature increases, modes with smaller 
characteristic lengths become unstable, and rearrangements can occur on 
shorter length scales. Hence this model predicts the presence of a growing 
length scale l(T) that converges toward I* when the temperature decreases 
toward T r . Growing dynamical length scales were observed numerically as 
the temperature decreases, see e.g. |116j and ref. therein. The curve l(T) 
could be measured by pinning the particles at the boundary of subsystems 
of size I, and by considering the dependence of the structural relaxation time 
with temperature for given I. Note that a similar test was already proposed 
in |119j to test dynamical length scales. 

When the temperature decreases below T r , we expect activated events 
to relax the glass structure. This is plausibely enhanced by the neighboring 
elastic instability, and by the large amount of nearly unstable modes. Finding 
the relaxation process of such weak structures is a necessary next step. Such 
models may lead to the prediction for the super-activated dependence of the 
relaxation time when the temperatures decreases and the elastic instability 
goes away. Once again, rigidity is not a local criterion, but demands the 
existence of a relation between coordination number and pressure on length 
scale I*. Hence it is not excluded that the fluctuations of quantities such as 
pressure and coordination on such distances may trigger relaxations. 
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10.2.3 Granular matter 

There are open questions on the way force propagates in amorphous sys- 
tems such as granular matter. At large distances, a granular pack should 
behave as a continuous elastic medium, at least in a small linear regime. As 
we discussed in Chapter EJ the situation is different in anisotropic isostatic 
systems, whose elastic behaviors can be described by hyperbolic equations. 
Nevertheless, there is as yet no description of force propagation in isotropic 
isostatic systems. More importantly, we do not know how force propagation 
evolves toward a normal continuous elastic behavior when the coordination 
number increases, or when friction is present. It was proposed in j^Z] that 
in anisotropic system, normal elasticity is recovered above /*, the distance 
at which soft modes disappear. Our derivation of the anomalous modes sug- 
gests that the transition toward continuous elasticity could occur at distances 
shorter than I*. In particular for a force dipole, or a local strain, as we dis- 
cussed in Chapter |H] it is plausible that the characteristic transverse length 
l t ~ 5z~2 characterizes the cross-over isostaticity/elasticity. Note that this 
does not preclude that the response to a force mono-pole has a different cross- 
over. It would be of much interest to investigate these subtle properties. It 
could be done in simulations of elastic spheres near the jamming threshold. 
Experimentally such tests would require to find efficient ways to modulate 
the coordination and the distance from isostaticity |118j . 

The length scales we are talking about obviously depend on the system 
considered, but we expect them to be typically of the order of a few tens 
of particle sizes. Thus they may not affect for example the building of sand 
castles, as a continuum description is expected to be valid for macroscopic 
objects. Nevertheless the properties of an assembly of grains at such length 
scales might play a crucial role in the rheology of these systems, both in the 
solid and in the dense liquid phases. Similar length scales were observed in 
the spatial correlations of dense granular flows on a slope, and are probably 
related to the surprising dependence of the thickness h of the flow with the 
angle of the sloped ^7|. Bending a layer of sand toward its avalanche angle 
is equivalent to imposing a shear. Thus understanding how shear affects 
the anomalous modes might shed light on this problem. The presence of 
a fixed boundary at a distance h of the free surface certainly increases the 
anomalous modes frequency. Making the most simple assumptions that (i) 
the shear decreases linearly the characteristic energy of the anomalous modes 
and (ii) the presence of of a fixed boundary is equivalent to an increase of 
coordination of order 1/h gives uj\ m ~ A\B(5z + 1/h) — A3 sin 9. Imposing 
uj am = yields a dependence of h of the form h{6) ~ jzg~: which is not 
inconsistent with the empirical data. 
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Another interesting property of granular matter is compaction: if a con- 
tainer of sand is tapped from below, the density of the system increases. The 
process displays many glassy features: the dynamics becomes very slow with 
time, aging and memory effect are observed [16J. When a granular pack is 
tapped weakly enough to avoid fluidization, there are two main causes for 
the irreversible events that lead to compaction. One the one hand, contacts 
on the Coulomb cone can slide. On the other hand, the destabilizing ef- 
fect of the pressure wave can cause the buckling of anomalous modes. This 
structural buckling increases the coordination and the packing fraction, as 
sketched in Fig. ()6.2|) . Hence the structure becomes more and more stable, 
and the compaction dynamics slows down. As the coordination rises the 
length scale I* decreases. This may be possible to test this presiction using 
X-ray microtomography. 
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